Version:

This documentation is not for the latest stable Salvus version.

This tutorial is presented as Python code running inside a Jupyter Notebook, the recommended way to use Salvus. To run it yourself you can copy/type each individual cell or directly download the full notebook, including all required files.

Data Adaptive Mesh Masking

In many seismologic applications, the source reveiver distribution does not follow a simple geometry such as a square or a circle. Due to the fact that sources are mostly confined to plate boundaries and stations are almost always on land, very complicated domain shapes are common. See [1] for an application.

In this tutorial, we show two different ways of adapting a mesh to the source receiver distribution my removing elements from a larger mesh that are not passed by waves of interest. The example we use here is a quake in Turkey recorded on the Search Results USArray Reference Network (_US-REF).

Copy
%matplotlib inline
%config Completer.use_jedi = False

from salvus.mesh.simple_mesh import SmoothieSEM, Globe3D
from salvus import namespace as sn

from obspy.clients.fdsn import Client

# a quake in Turkey that we will use in this tutorial, original data from IRIS spud:
# http://service.iris.edu/fdsnws/event/1/query?eventid=2847365
# http://ds.iris.edu/spudservice/momenttensor/gcmtid/C201003241411A/quakeml#momenttensor
source = sn.simple_config.source.seismology.SideSetMomentTensorPoint3D(
    latitude=38.82,
    longitude=40.14,
    depth_in_m=4500,
    side_set_name="r1",
    mrr=5.47e15,
    mtt=-4.11e16,
    mpp=3.56e16,
    mrt=2.26e16,
    mrp=-2.25e16,
    mtp=1.92e16,
)

# get USarray stations from iris
inv = Client("IRIS").get_stations(
    network="_US-REF", level="station", format="text"
)

receivers = sn.simple_config.receiver.seismology.parse(
    inv, dimensions=3, fields=["displacement"]
)

# prepare an event collection that will later be used to mask the mesh to a region of interest
event_collection = sn.EventCollection.from_sources(
    sources=[source], receivers=receivers
)

Method 1: Surface Based

In the first approach, we only mask out element based on there lateral position on the sphere and ignore the depth (this approach is used in [1]). To achieve this, a covex hull is built from the sources and receivers and all elements within that hull as well as those within a specified distance are retained in the mesh.

from salvus.mesh.mask_generators import SurfaceMaskGenerator

sm = Globe3D()
sm.basic.model = "prem_iso_one_crust"
sm.basic.min_period_in_seconds = 100.0
sm.basic.elements_per_wavelength = 2.0
sm.spherical.min_radius = 4000.0


# use event collection to create a surface mask
smg = SurfaceMaskGenerator(
    event_collection,
    number_of_points=1000,
    distance_in_km=1000.0,
)

# hand over the mask as a callback funcion
sm.create_mesh(mesh_processing_callback=smg)
/miniconda/envs/salvus/lib/python3.7/_collections_abc.py:702: MatplotlibDeprecationWarning: The global colormaps dictionary is no longer considered public API.
  return len(self._mapping)
/miniconda/envs/salvus/lib/python3.7/_collections_abc.py:720: MatplotlibDeprecationWarning: The global colormaps dictionary is no longer considered public API.
  yield from self._mapping
When run interactively in a Jupyter Notebook, this output cell contains a widget. Click the button below to load a preview of it. Please note that most interactive functionality does not work without a running Python kernel.
<salvus.mesh.unstructured_mesh.UnstructuredMesh at 0x7f3b57386190>

Alternatively, we can use rays computed with obspy.TauPy to mask out all elements beyond a specified distance from the rays. E.g. for teleseismic body wave tomography, this reduces the number of elements further in comparison to the surface based masking. This may be particularly interesting in the context of SmoothieSEM meshes, however, this works identically with standard cubed sphere type meshes.

In this particular example we only use the P-waves, as would be applicable in a teleseismic P-wave tomography:

from salvus.mesh.mask_generators import RayMaskGenerator

sm = SmoothieSEM()
sm.basic.model = "prem_iso_one_crust"
sm.basic.min_period_in_seconds = 50.0
sm.basic.elements_per_wavelength = 2.0
sm.basic.number_of_lateral_elements = 12
sm.spherical.min_radius = 1500.0

# no refinement, rotation to a source in Turkey
sm.basic.number_of_lateral_elements = 10
sm.source.latitude = source._initial_arguments["latitude"]
sm.source.longitude = source._initial_arguments["longitude"]

# use event collection to create a ray mask
rmg = RayMaskGenerator(
    event_collection,
    phases=["P"],
    number_of_points_per_ray=100,
    distance_in_km=1000.0,
)

sm.create_mesh(mesh_processing_callback=rmg)
/miniconda/envs/salvus/lib/python3.7/_collections_abc.py:702: MatplotlibDeprecationWarning: The global colormaps dictionary is no longer considered public API.
  return len(self._mapping)
/miniconda/envs/salvus/lib/python3.7/_collections_abc.py:720: MatplotlibDeprecationWarning: The global colormaps dictionary is no longer considered public API.
  yield from self._mapping
When run interactively in a Jupyter Notebook, this output cell contains a widget. Click the button below to load a preview of it. Please note that most interactive functionality does not work without a running Python kernel.
<salvus.mesh.unstructured_mesh.UnstructuredMesh at 0x7f3b5d49e090>

[1] van Herwaarden, Dirk Philip, Christian Boehm, Michael Afanasiev, Solvi Thrastarson, Lion Krischer, Jeannot Trampert, and Andreas Fichtner. 2020. “Accelerated Full-Waveform Inversion Using Dynamic Mini-Batches.” Geophysical Journal International 221 (2): 1427–38. https://doi.org/10.1093/gji/ggaa079.

PAGE CONTENTS