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.

# Marmousi Model

## Introduction

The AGL Elastic Marmousi model is the most recent generation of the classic Marmousi seismic benchmark. In this tutorial we will cover how one can use this model within Salvus. In particular, we will cover:
• How to read a SEG-Y file using Obspy
• How to read in a rectilinear model using xarray
• How to interpolate a 2-D regularly-gridded model onto the mesh
• How to generate a synthetic seismic shotgather using the model

## Environment setup

As always, our first job is to import the relevant Python packages into our notebook.
Copy
%matplotlib inline
%config Completer.use_jedi = False

from pathlib import Path
import os
import requests
import shutil
from typing import List

import scipy
import obspy
import numpy as np
import xarray as xr
import matplotlib.pyplot as plt

from salvus.flow import api
from salvus.flow import simple_config as config

from salvus.toolbox import toolbox

SALVUS_FLOW_SITE_NAME = os.environ.get("SITE_NAME", "local")

## Model building

Information about the Marmousi model we'll be using can be found here. The model is stored as a SEG-Y file, and the next cell will download and unpack the model for you. Note that the model file is relatively large (~150 MB), so make sure you have enough disk space before continuing.
# Download Marmousi model ~150 MB.
target_path = "elastic-marmousi-model.tar.gz"
if not Path(target_path).exists():
url = (
"https://s3.amazonaws.com/open.source.geoscience/open_data"
"/elastic-marmousi/elastic-marmousi-model.tar.gz"
)

response = requests.get(url, stream=True)
if response.status_code == 200:
with open(target_path, "wb") as f:
shutil.unpack_archive(target_path)

### Reading the model

With our environment is set up, and the model downloaded, we can begin building the Salvus mesh. Most material models in seismic exploration are discretized onto a regular grid, and the open-source SalvusPythonToolkit provides a simple interface to create two or three dimensional Salvus meshes based on such models. As models can be stored in a wide variety of format such as SEG-Y, or RSF, we use xarray as an intermediary to encapsulate both the model parameters and geometry. More info on xarray, including extensive documentation, can be found here. Reading our Marmousi model into an xarray.Dataset is trivial, and the process can be inspected by opening the acompanying getting_started_tools.py file.
In general, if you are working working with a regularly gridded model in either two or three dimensions, a recipie similar to the function provided should be all you need to do to get your model into a form Salvus Mesh can understand. All that is required is an xarray dataset with the following information
• The coordinate values in meters (named 'x', 'y', and optionally 'z')
• The material parameters in SI units, with names following the parameter naming convention
# Read marmousi model.
model_directory = Path("./elastic-marmousi-model/model")
marmousi_model = toolbox.read_elastic_marmousi(model_directory, ["VP", "RHO"])
A nice feature of xarray.Dataset object is the ability to quickly plot the contained data. Let's do this below as a way to quickly QC that the model reading went correctly.
marmousi_model["VP"].T.plot(aspect="auto", figsize=(15, 5))
plt.show()
marmousi_model["RHO"].T.plot(aspect="auto", figsize=(15, 5))
plt.show()
Yup, all looks good!

### Meshing parameters

#### Frequency content

When building material models for use with wave propagation solvers such as Salvus, there are a few factors one needs to consider. Obviously we need to ensure that the spatial complexity of the model is adequetly captured by whatever spatial discretization we are choosing. What is less obvious however is that the details of the solver's spatial discretization are intimitely to both the spatial complexity of the model and the spatial complexity of the physical process we are solving for. When considering the propagation of seismic waves we are somewhat lucky that the solutions are band-limited and the spatial complexity of the solution can be reliably estimated before any simulations are run.
Using the simple formula $\lambda_{\text{min}} = v_{\text{min}} \cdot f_{\text{max}}^{-1}$, we have an estimate $\lambda_{\text{min}}$ for the minimum wavelength in our simulation when given the minimum velocity in the model ($v_{\text{min}}$