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.

Surface topography

The effects of free-surface topography on seismic waveforms are significant, and any local study considering higher-frequency data must take topography into account. In this tutorial we'll explore how to add realistic free-surface topography to a SimulationConfiguration object, and investigate how high-order accurate topography / material interpolation can be used to improve interpolation performance. To accomplish this, we'll focus on a real-life use case using the area around Mt. St. Helens in Washington state as an example.

As always, we'll start by importing the Salvus namespace.

Copy
%matplotlib inline
import os
import pathlib
import salvus.namespace as sn
# This notebook will use this variable to determine which
# remote site to run on.
SALVUS_FLOW_SITE_NAME = os.environ.get("SITE_NAME", "local")
PROJECT_DIR = "project"

Getting topography files

Salvus natively supports surface topography elevation models available via the AppEEARS webservice. The AppEEARS project is run by the USGS, and offers "a simple and efficient way to access and transform geospatial data from a variety of federal data archives". For our purposes we're interested in data stored from NASA's Shuttle Radar Topography Mission, or SRTM. This project provides elevation measurements with reference to the WGS84 ellipsoid at 30m resolution for most of the globe. Access to the data is provided through a simple graphical user interface, as shown in the figure below -- it is this very request which we will use in the tutorial today. To use AppEEARS yourself you simply need to create an account at the website and get going! The workflow we present below is generic, and can be applied to any UTM domain for which SRTM data is available.

Once an AppEEARS request is made, the data will be fetched and you will get an email informing you that it is ready for download. The download composes of at least two files:

  • A .json file containing the parameters of your request, as well as the spatial domain for which it was made
  • A NetCDF file containing the SRTM data at the resolution you requested.

The Domain constructor within Salvus project can recognize AppEEARS request files, and can use them to set automatically set up a cartesian domain in the correct UTM zone. We do this below. You'll notice that we supply an extra parameter: shrink_domain. This factor exists to make it easy to isotropically shrink the domain specified in the AppEEARS request. Specifying such a value is necessary in our case as the mesh we create will eventually require absorbing boundary layers to be attached. If these layers extend past the boundaries of the topography files, we should expect some artefacts in the interpolation of the topography. Shrinking the domain here provides an easy way to ensure that the topography model is still valid for the absorbing region.

d = sn.domain.dim3.UtmDomain.from_appeears_request(
    json_file="./appeears.json", shrink_domain=10000.0
)

Once our domain is created, we can initialize our project as before.

# Uncomment the following line to delete a
# potentially existing project for a fresh start
# !rm -rf proj
if pathlib.Path(PROJECT_DIR).exists():
    print("Opening existing project.")
    p = sn.Project(path="project")
else:
    print("Creating new project.")
    p = sn.Project.from_domain(path=PROJECT_DIR, domain=d)
Creating new project.

Just as the domain has a constructor which can take an AppEEARS request, the actual topography data can also be initalized using a similar interface. Here we use the second file passed to us from AppEEARS, which is the NetCDF file containing the actual topography data. This inteface also takes a few additional parameters which determine how the topography is read and constructed. Here we pass decimate_topo_factor which only reads every "nth" value in the data file. Since we will be working with relatively low-frequency simulations we will not need topography data every 30m. Decimating the data allows us to work with a smaller file in memory with almost no effect on the final mesh. We also pass the utm object stored within the Domain. The object is a pyproj.Proj projection object which defines how we transform the coordinates in the APPEEARS request to those in our UTM domain.

st = sn.topography.cartesian.SurfaceTopography.from_appeears_request(
    name="topo", data="./sthelens.nc", decimate_topo_factor=100, utm=d.utm
)

The interface we just used above to read in topography data is very general, and we'll see in a future tutorial how to read in surface, Moho, or discontinuity topography using a similar simplified interface. As is true for the volumetric models used in Salvus, once the topography files are read and processed they are stored as xarray datasets in memory. These datasets are simple to manipulate and visualize, and we can for instance plot the digital elevation model in the notebook by simply plotting the xarray dataset.

st.ds.dem.plot()
<matplotlib.collections.QuadMesh at 0x7fc7893abb90>

Looks like St. Helens to me! As with all data we want to keep, we'll add the processed topography file to the project in the standard way.

p += st

Now we'll add an explosive source and a line of receivers to the project. Here we use the SideSet source and receiver classes. The functionality of these classes is documented here. Essentially, they help when the actual surface of the domain is non-trivial, as is the case in this example. When these sources and receivers are associated with a mesh, a small optimization problem is solved to ensure that they are placed exactly at the surface of the deformed mesh, plus or minus any offset one might specify. Here we specify the entities in UTM coordinates, place the source 1000 meters below the free surface, and place a line of 50 receivers directly at the free surface.

# Src / Rec reference coordinates.
src_x, src_y, src_z = 562700.0, 5112500.0, 4000.0
rec_x, rec_y, rec_z = 564700.0, 5115500.0, 4000.0
# Place an explosive source 1000 m below the free surface.
src = sn.simple_config.source.cartesian.SideSetMomentTensorPoint3D(
    point=(src_x, src_y, src_z),
    direction=(0, 0, 1),
    side_set_name="z1",
    mxx=1e21,
    myy=1e21,
    mzz=1e21,
    myz=0.0,
    mxz=0.0,
    mxy=0.0,
    offset=-1000.0,
)
# Place a line of receivers at the free surface.
rec = [
    sn.simple_config.receiver.cartesian.SideSetPoint3D(
        point=(rec_x, rec_y + _i * 100, rec_z),
        direction=(0, 0, 1),
        side_set_name="z1",
        fields=["velocity"],
        station_code=f"XX_{_i}",
    )
    for _i in range(50)
]
# Add the event collection to the project.
p += sn.EventCollection.from_sources(sources=[src], receivers=rec)

The next few commands should now be familiar after the previous tutorials. Below we initialize an isotropic elastic background model, a homogeneous model configuration constructed from this background model and our event configuration. We could also consider adding a 1-D background model, or a 3-D material model, at this stage, but we'll save this for a future tutorial.

bm = sn.model.background.homogeneous.IsotropicElastic(
    rho=2600.0, vp=3000.0, vs=1875.5
)
mc = sn.model.ModelConfiguration(background_model=bm)
ec = sn.EventConfiguration(
    wavelet=sn.simple_config.stf.Ricker(center_frequency=1.0),
    waveform_simulation_configuration=sn.WaveformSimulationConfiguration(
        end_time_in_seconds=6.0
    ),
)

Below we do something new, but which follows the model API we're now familiar with. In the cell below we initialize a TopographyConfiguration object. As is true with the ModelConfiguration, this object takes as arguments one or more topography models. A background topography model does not have much of a meaning, so it is not defined for topography configurations.

tc = sn.topography.TopographyConfiguration(topography_models="topo")

Next we initialize an object which will inform the extrusion of the domain for the absorbing boundaries.

ab = sn.AbsorbingBoundaryParameters(
    reference_velocity=3000.0,
    reference_frequency=5.0,
    number_of_wavelengths=3.5,
)

And now we've got enough information to define our complete simulation configuration. We create the configuration below and add it to the project as usual.

p += sn.SimulationConfiguration(
    name="sim_1st_order_topo",
    tensor_order=1,
    model_configuration=mc,
    event_configuration=ec,
    absorbing_boundaries=ab,
    elements_per_wavelength=1,
    max_depth_in_meters=2000.0,
    max_frequency_in_hertz=2.0,
    topography_configuration=tc,
)

We can now visualize the simulation mesh with our topography on top. Looks pretty fancy!

p.viz.nb.simulation_setup("sim_1st_order_topo", ["event_0000"])
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.
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.flow.simple_config.simulation.Waveform object at 0x7fc76ced0b10>

It's now time to run a simulation and analyze the results. Since this is a larger 3-D simulation it will not run as quickly as the Lamb's problem tutorial, so you may have to wait a moment.

p.simulations.launch(
    "sim_1st_order_topo",
    events=["event_0000"],
    site_name=SALVUS_FLOW_SITE_NAME,
    ranks_per_job=2,
    max_block_in_seconds=240,
)
Uploading 1 files...

🚀  Submitted [email protected]local
After 0.1 seconds: running. Sleeping for 0.2500 seconds before checking job status again.
After 0.4 seconds: running. Sleeping for 0.2500 seconds before checking job status again.
After 0.7 seconds: running. Sleeping for 0.2500 seconds before checking job status again.
After 0.9 seconds: running. Sleeping for 0.2500 seconds before checking job status again.
After 1.2 seconds: running. Sleeping for 0.2500 seconds before checking job status again.
After 1.5 seconds: running. Sleeping for 0.2500 seconds before checking job status again.
After 1.8 seconds: running. Sleeping for 0.2500 seconds before checking job status again.
After 2.0 seconds: running. Sleeping for 0.2500 seconds before checking job status again.
After 2.3 seconds: running. Sleeping for 0.2500 seconds before checking job status again.
After 2.6 seconds: running. Sleeping for 0.2500 seconds before checking job status again.
After 2.8 seconds: running. Sleeping for 0.2500 seconds before checking job status again.
After 3.1 seconds: running. Sleeping for 0.2500 seconds before checking job status again.
After 3.3 seconds: running. Sleeping for 0.2500 seconds before checking job status again.
After 3.6 seconds: running. Sleeping for 0.2500 seconds before checking job status again.
After 3.8 seconds: running. Sleeping for 0.2500 seconds before checking job status again.
After 4.1 seconds: running. Sleeping for 0.2500 seconds before checking job status again.
After 4.4 seconds: running. Sleeping for 0.2500 seconds before checking job status again.
After 4.6 seconds: running. Sleeping for 0.2500 seconds before checking job status again.
After 4.9 seconds: running. Sleeping for 0.2500 seconds before checking job status again.
After 5.1 seconds: running. Sleeping for 0.2500 seconds before checking job status again.
After 5.4 seconds: running. Sleeping for 0.2500 seconds before checking job status again.
After 5.6 seconds: running. Sleeping for 0.2500 seconds before checking job status again.
After 5.9 seconds: running. Sleeping for 0.2500 seconds before checking job status again.
After 6.1 seconds: running. Sleeping for 0.2500 seconds before checking job status again.
After 6.4 seconds: running. Sleeping for 0.2500 seconds before checking job status again.
After 6.7 seconds: running. Sleeping for 0.2500 seconds before checking job status again.
After 6.9 seconds: running. Sleeping for 0.2500 seconds before checking job status again.
After 7.2 seconds: running. Sleeping for 0.2500 seconds before checking job status again.
After 7.5 seconds: running. Sleeping for 0.2500 seconds before checking job status again.
After 7.7 seconds: running. Sleeping for 0.2500 seconds before checking job status again.
After 8.0 seconds: running. Sleeping for 0.2500 seconds before checking job status again.
After 8.2 seconds: running. Sleeping for 0.2500 seconds before checking job status again.
After 8.6 seconds: running. Sleeping for 0.2500 seconds before checking job status again.
After 8.8 seconds: running. Sleeping for 0.2500 seconds before checking job status again.
After 9.1 seconds: running. Sleeping for 0.2500 seconds before checking job status again.
After 9.3 seconds: running. Sleeping for 0.2500 seconds before checking job status again.
After 9.6 seconds: running. Sleeping for 0.2500 seconds before checking job status again.
After 9.9 seconds: running. Sleeping for 0.2500 seconds before checking job status again.
After 10.1 seconds: running. Sleeping for 0.2500 seconds before checking job status again.
After 10.5 seconds: running. Sleeping for 0.2500 seconds before checking job status again.
After 10.7 seconds: running. Sleeping for 0.2500 seconds before checking job status again.
After 11.0 seconds: running. Sleeping for 0.2500 seconds before checking job status again.
After 11.2 seconds: running. Sleeping for 0.2500 seconds before checking job status again.
After 11.5 seconds: running. Sleeping for 0.2500 seconds before checking job status again.
After 11.8 seconds: running. Sleeping for 0.2500 seconds before checking job status again.
After 12.1 seconds: running. Sleeping for 0.2500 seconds before checking job status again.
After 12.3 seconds: running. Sleeping for 0.2500 seconds before checking job status again.
After 12.6 seconds: running. Sleeping for 0.2500 seconds before checking job status again.
After 12.9 seconds: running. Sleeping for 0.2500 seconds before checking job status again.
After 13.1 seconds: running. Sleeping for 0.2500 seconds before checking job status again.
After 13.4 seconds: running. Sleeping for 0.2500 seconds before checking job status again.
After 13.7 seconds: running. Sleeping for 0.2500 seconds before checking job status again.
After 14.0 seconds: running. Sleeping for 0.2500 seconds before checking job status again.
After 14.2 seconds: running. Sleeping for 0.2500 seconds before checking job status again.
After 14.5 seconds: running. Sleeping for 0.2500 seconds before checking job status again.
After 14.8 seconds: running. Sleeping for 0.2500 seconds before checking job status again.
After 15.1 seconds: running. Sleeping for 0.2500 seconds before checking job status again.
After 15.5 seconds: running. Sleeping for 0.2500 seconds before checking job status again.
After 15.8 seconds: running. Sleeping for 0.2500 seconds before checking job status again.
1
p.simulations.query()
True

And as before when the simulation is complete we can fetch the data and create a shotgather. Given that we have a relatively dense sampling of receivers, looking at both the wiggles and the shotgather together can be nice.

p.waveforms.get(data_name="sim_1st_order_topo", events=["event_0000"])[0].plot(
    component="X",
    receiver_field="velocity",
    plot_types=["wiggles", "shotgather"],
)

We of course want to do something a bit more than just run a simulation here -- we're going to breifly investigate the effect of topography on our waveforms. If you remember from the theoretical background, Salvus allows one to use higher-order elements to discretize both the model and topography. Below we create a new simulation configuration where the only difference is that we increase the order of topography interpolation to 4. Let's run this and see the results.

# Create the configuration.
p += sn.SimulationConfiguration(
    name="sim_4th_order_topo",
    tensor_order=4,
    model_configuration=mc,
    event_configuration=ec,
    absorbing_boundaries=ab,
    elements_per_wavelength=1,
    max_depth_in_meters=2000.0,
    max_frequency_in_hertz=2.0,
    topography_configuration=tc,
)

# Launch the simulation.
p.simulations.launch(
    "sim_4th_order_topo",
    events=["event_0000"],
    site_name=SALVUS_FLOW_SITE_NAME,
    ranks_per_job=2,
    max_block_in_seconds=180,
)

# Visualize
p.viz.nb.simulation_setup("sim_4th_order_topo", ["event_0000"])
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.
Uploading 1 files...

🚀  Submitted [email protected]local
After 0.1 seconds: running. Sleeping for 0.2500 seconds before checking job status again.
After 0.4 seconds: running. Sleeping for 0.2500 seconds before checking job status again.
After 0.7 seconds: running. Sleeping for 0.2500 seconds before checking job status again.
After 1.0 seconds: running. Sleeping for 0.2500 seconds before checking job status again.
After 1.3 seconds: running. Sleeping for 0.2500 seconds before checking job status again.
After 1.6 seconds: running. Sleeping for 0.2500 seconds before checking job status again.
After 1.8 seconds: running. Sleeping for 0.2500 seconds before checking job status again.
After 2.1 seconds: running. Sleeping for 0.2500 seconds before checking job status again.
After 2.4 seconds: running. Sleeping for 0.2500 seconds before checking job status again.
After 2.7 seconds: running. Sleeping for 0.2500 seconds before checking job status again.
After 3.0 seconds: running. Sleeping for 0.2500 seconds before checking job status again.
After 3.3 seconds: running. Sleeping for 0.2500 seconds before checking job status again.
After 3.5 seconds: running. Sleeping for 0.2500 seconds before checking job status again.
After 3.8 seconds: running. Sleeping for 0.2500 seconds before checking job status again.
After 4.0 seconds: running. Sleeping for 0.2500 seconds before checking job status again.
After 4.4 seconds: running. Sleeping for 0.2500 seconds before checking job status again.
After 4.6 seconds: running. Sleeping for 0.2500 seconds before checking job status again.
After 4.9 seconds: running. Sleeping for 0.2500 seconds before checking job status again.
After 5.1 seconds: running. Sleeping for 0.2500 seconds before checking job status again.
After 6.0 seconds: running. Sleeping for 0.2500 seconds before checking job status again.
After 6.3 seconds: running. Sleeping for 0.2500 seconds before checking job status again.
After 6.5 seconds: running. Sleeping for 0.2500 seconds before checking job status again.
After 6.8 seconds: running. Sleeping for 0.2500 seconds before checking job status again.
After 7.0 seconds: running. Sleeping for 0.2500 seconds before checking job status again.
After 7.3 seconds: running. Sleeping for 0.2500 seconds before checking job status again.
After 7.5 seconds: running. Sleeping for 0.2500 seconds before checking job status again.
After 7.8 seconds: running. Sleeping for 0.2500 seconds before checking job status again.
After 8.1 seconds: running. Sleeping for 0.2500 seconds before checking job status again.
After 8.4 seconds: running. Sleeping for 0.2500 seconds before checking job status again.
After 8.6 seconds: running. Sleeping for 0.2500 seconds before checking job status again.
After 8.9 seconds: running. Sleeping for 0.2500 seconds before checking job status again.
After 9.1 seconds: running. Sleeping for 0.2500 seconds before checking job status again.
After 9.4 seconds: running. Sleeping for 0.2500 seconds before checking job status again.
After 9.6 seconds: running. Sleeping for 0.2500 seconds before checking job status again.
After 9.9 seconds: running. Sleeping for 0.2500 seconds before checking job status again.
After 10.2 seconds: running. Sleeping for 0.2500 seconds before checking job status again.
After 10.5 seconds: running. Sleeping for 0.2500 seconds before checking job status again.
After 10.7 seconds: running. Sleeping for 0.2500 seconds before checking job status again.
After 11.0 seconds: running. Sleeping for 0.2500 seconds before checking job status again.
After 11.3 seconds: running. Sleeping for 0.2500 seconds before checking job status again.
After 11.5 seconds: running. Sleeping for 0.2500 seconds before checking job status again.
After 11.8 seconds: running. Sleeping for 0.2500 seconds before checking job status again.
After 12.0 seconds: running. Sleeping for 0.2500 seconds before checking job status again.
After 12.4 seconds: running. Sleeping for 0.2500 seconds before checking job status again.
After 12.6 seconds: running. Sleeping for 0.2500 seconds before checking job status again.
After 12.9 seconds: running. Sleeping for 0.2500 seconds before checking job status again.
After 13.1 seconds: running. Sleeping for 0.2500 seconds before checking job status again.
After 13.4 seconds: running. Sleeping for 0.2500 seconds before checking job status again.
After 13.6 seconds: running. Sleeping for 0.2500 seconds before checking job status again.
After 13.9 seconds: running. Sleeping for 0.2500 seconds before checking job status again.
After 14.2 seconds: running. Sleeping for 0.2500 seconds before checking job status again.
After 14.5 seconds: running. Sleeping for 0.2500 seconds before checking job status again.
After 14.7 seconds: running. Sleeping for 0.2500 seconds before checking job status again.
After 15.0 seconds: running. Sleeping for 0.2500 seconds before checking job status again.
After 15.3 seconds: running. Sleeping for 0.2500 seconds before checking job status again.
After 15.5 seconds: running. Sleeping for 0.2500 seconds before checking job status again.
After 15.8 seconds: running. Sleeping for 0.2500 seconds before checking job status again.
After 16.0 seconds: running. Sleeping for 0.2500 seconds before checking job status again.
After 16.3 seconds: running. Sleeping for 0.2500 seconds before checking job status again.
After 16.6 seconds: running. Sleeping for 0.2500 seconds before checking job status again.
After 16.8 seconds: running. Sleeping for 0.2500 seconds before checking job status again.
After 17.1 seconds: running. Sleeping for 0.2500 seconds before checking job status again.
After 17.4 seconds: running. Sleeping for 0.2500 seconds before checking job status again.
After 17.6 seconds: running. Sleeping for 0.2500 seconds before checking job status again.
After 17.9 seconds: running. Sleeping for 0.2500 seconds before checking job status again.
After 18.1 seconds: running. Sleeping for 0.2500 seconds before checking job status again.
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.flow.simple_config.simulation.Waveform object at 0x7fc76d14f610>

Now we can use the comparison GUI to see how the more accurate interpolation of topography affected our waveforms. Quite a bit as you can see!

p.simulations.query()
p.viz.nb.waveforms(
    ["sim_1st_order_topo", "sim_4th_order_topo"], receiver_field="velocity"
)
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.
# Create the configuration.
# p += sn.SimulationConfiguration(
#     name="sim_1st_order_topo_large",
#     tensor_order=1,
#     model_configuration=mc,
#     event_configuration=ec,
#     absorbing_boundaries=ab,
#     elements_per_wavelength=1,
#     max_depth_in_meters=2000.0,
#     max_frequency_in_hertz=5.0,
#     topography_configuration=tc,
# )
# if False:
#     p.simulations.launch(
#         "sim_1st_order_topo_large",
#         events=["event_0000"],
#         site_name="cluster",
#         ranks_per_job=36,
#         wall_time_in_seconds_per_job=600,
#     )
# if False:
#     p.waveforms.get(
#         data_name="sim_1st_order_topo_large", events=["event_0000"]
#     )[0].plot(["wiggles", "shotgather"])
PAGE CONTENTS