Version:

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 by 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
from salvus import namespace as sn

# 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,
)

receivers = sn.simple_config.receiver.seismology.parse(
    "data/inventory.xml", 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 = sn.simple_mesh.basic_mesh.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(processing_function=smg)
<salvus.mesh.data_structures.unstructured_mesh.unstructured_mesh.UnstructuredMesh object at 0x7f506910ed90>
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.
In this particular example we only use the P-waves, as would be applicable in a teleseismic P-wave tomography:
from salvus.mesh.algorithms.mask_generators import RayMaskGenerator

sm = sn.simple_mesh.basic_mesh.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 = 1500.0

# 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(processing_function=rmg)
<salvus.mesh.data_structures.unstructured_mesh.unstructured_mesh.UnstructuredMesh object at 0x7f50572657d0>
[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