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.

Global and Regional Seismic Wave Propagation with Smothiesem

This tutorial demonstrates how to prepare meshes for global seismic wave simulations using the anisotropic adaptive mesh refinement rechnique (aamr), which in this context is also known as smoothiesem [1, 2]. The meshes are specifically build for a particular source location to take advantage of the lateral smoothness of the wavefield and reduce the numerical burden significantly for global and regional scale simulations.

The mesh can be refined in predefined depth, distance and azimuth regions. Additionally, the domain can be restricted to a region of interest as discussed in the separate data adaptive mesh masking tutorial.

In this tutorial we demonstrate the basic concept and most important parameters to build these meshes.

Copy
%matplotlib inline
%config Completer.use_jedi = False

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

from obspy.clients.fdsn import Client

period_global = 100.0

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

Full Sphere Meshes

First, we build some full sphere meshes to highlight the different lateral refinement options.

1) no lateral refinement

The simplemost smoothiesem mesh is a full sphere with no lateral refinements and the symmetry axis is aligned with the seismic source location. The tensor order controls boths the accuracy of representing the seismic velocities as well as the geometry of the domain (here: the sphere). On top, the parameter nlat allows to vary the number of elements in the lateral direction that is used throughout the whole mesh. Compare Figure 9 in [1] to choose values appropriate for your application.

sm = SmoothieSEM()
sm.basic.model = "prem_iso_one_crust"
sm.basic.min_period_in_seconds = period_global
sm.basic.elements_per_wavelength = 2.0

sm.advanced.tensor_order = 2
sm.basic.number_of_lateral_elements = 4
sm.source.latitude = source._initial_arguments["latitude"]
sm.source.longitude = source._initial_arguments["longitude"]

sm.create_mesh()
/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 0x7fbf11d20d10>

This type of refinement is particularly useful to increase resolution of the surface towards the equator, hence to improve resolution of surface waves as well as surface topography. Here, theta is the angular distance from the source in the range [0°, 180°]. Multiple such refinements can be combined, but refinement boundaries should not cross to ensure high quality elements: each refinement should be fully contained in all previous refinements.

sm.refinement.lateral_refinements = [
    {
        "theta_min": 60.0,
        "theta_max": 150.0,
        "r_min": 6000.0,
    },
    {
        "theta_min": 80.0,
        "theta_max": 130.0,
        "r_min": 6200.0,
    },
]

sm.create_mesh()
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 0x7fbf12ae6990>

This type of refinement may be used to increase resolution in particular areas of interest defined by a distance and azimuth range [-180°, 180°] for [phi_min, phi_max] relative to the source azimuths measured clockwise from north. This may be a practical alternative to absorbing boundaries and a azimuthally constrained domain.

sm.refinement.lateral_refinements = [
    {
        "theta_min": 40.0,
        "theta_max": 110.0,
        "r_min": 6000.0,
        "phi_min": -45.0,
        "phi_max": 45.0,
    }
]

sm.source.azimuth = 270.0

sm.create_mesh()
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 0x7fc04ce42a10>

The lateral refinements can also be contrained to a depth region not at the surface to resolve e.g. small scale structure at a certain dept of interest. This can be achieved with the r_max parameter of ther lateral refinement object:

sm.refinement.lateral_refinements = [
    {
        "theta_min": 40.0,
        "theta_max": 110.0,
        "r_min": 3000.0,
        "r_max": 4000.0,
        "phi_min": -45.0,
        "phi_max": 45.0,
    }
]

sm.create_mesh()
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 0x7fbf12779c90>

The meshing region can be constrained in depth and distance to create regional or local meshes with higher resolution. This is provided by the following options in the simple mesh object. These parameters can freely be combined with the refinement parameters from above.

sml = SmoothieSEM()
sml.basic.model = "prem_iso_one_crust"
sml.basic.min_period_in_seconds = 10.0
sml.basic.elements_per_wavelength = 2.0

sml.spherical.min_radius = 6000.0
sml.chunk.max_colatitude = 10.0

sml.advanced.tensor_order = 2
sml.basic.number_of_lateral_elements = 4
sml.source.latitude = source._initial_arguments["latitude"]
sml.source.longitude = source._initial_arguments["longitude"]

sml.create_mesh()
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 0x7fbf13a0c890>

SmoothieSEM meshes are also capable to include smooth 3D stucture both in terms of the material properties (i.e. seismic velocities) as well as the domain shape (e.g. topography). However, special care needs to be taken to choose appropriate resolution for the mesh as well as the 3D structure and this is currently not automatized. Instead, it is the users responsibility to explicitly choose the resolution of the input data used in the meshing and verify not only it is represented without aliasing on the mesh, but also that the smoothieSEM approximation is still apropriate for the given application e.g. by a convergence test refining in azimuthal direction.

First, download the data files:

import pathlib
import requests

for fname in [
    "topography_earth2014_egm2008_lmax_256.nc",
    "bathymetry_earth2014_lmax_256.nc",
    "moho_topography_crust_1_0_egm2008.nc",
]:

    topography_file = pathlib.Path(fname)

    if not topography_file.exists():
        r = requests.get(
            f"https://data.mondaic.com/topography-data/{fname}",
            stream=True,
        )

        assert r.ok

        with topography_file.open("wb") as f:
            f.write(r.raw.read())

Topography, crustal thickness, ocean loading and ellipticity can be included in the meshing using the following options:

sm = SmoothieSEM()
sm.basic.model = "prem_iso_one_crust"
sm.basic.min_period_in_seconds = 100.0
sm.basic.elements_per_wavelength = 2.0
sm.basic.number_of_lateral_elements = 6
sm.advanced.tensor_order = 2

sm.source.latitude = 38.82
sm.source.longitude = 40.14

sm.topography.topography_file = "topography_earth2014_egm2008_lmax_256.nc"
sm.topography.topography_varname = (
    "topography_earth2014_egm2008_lmax_256_lmax_32"
)

sm.topography.moho_topography_file = "moho_topography_crust_1_0_egm2008.nc"
sm.topography.moho_topography_varname = (
    "moho_topography_crust_1_0_egm2008_lmax_32"
)

The topography files can be directly opened with paraview (using the NetCDF Reader) to vizualize the available datasets. The variable name is composed of the filename appended with the maximum resolution in terms of the spherical harmonic degree and order lmax. For reasonable accuracy, lmax should not be larger than approximately number_of_lateral_elements * tensor_order, so the settings in the example here are very optimistic to keep it small enough.

# WGS84 ellipticity value
sm.spherical.ellipticity = 0.0033528106647474805

m = sm.create_mesh()

# for vizualisation compute the radius (note that this in case includes the ellipticity)
m.attach_field("r", (m.points ** 2).sum(axis=1) ** 0.5)

# choose 'OCEANLOAD' or 'r' in the widget to see the 3D data
m
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 0x7fbf1823b750>

[1] Driel, Martin van, C Boehm, Lion Krischer, and Michael Afanasiev. 2020. “Accelerating Numerical Wave Propagation Using Wavefield Adapted Meshes. Part I: Forward and Adjoint Modelling.” Geophysical Journal International 221 (3): 1580–90. https://doi.org/10.1093/gji/ggaa058.

[2] Thrastarson, Solvi, Martin van Driel, Lion Krischer, Christian Boehm, Michael Afanasiev, Dirk-Philip van Herwaarden, and Andreas Fichtner. 2020. “Accelerating Numerical Wave Propagation by Wavefield Adapted Meshes. Part II: Full-Waveform Inversion.” Geophysical Journal International 221 (3): 1591–1604. https://doi.org/10.1093/gji/ggaa065.

PAGE CONTENTS