Version:

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
As always, our first job is to import the relevant Python packages into our notebook.
Copy
%matplotlib inline
%config Completer.use_jedi = False

import pathlib
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")
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 pathlib.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:
            f.write(response.raw.read())
        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 = pathlib.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!

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 λmin=vminfmax1\lambda_{\text{min}} = v_{\text{min}} \cdot f_{\text{max}}^{-1}, we have an estimate λmin\lambda_{\text{min}} for the minimum wavelength in our simulation when given the minimum velocity in the model (vminv_{\text{min}}) and the maximum frequency in the forcing function (fmaxf_{\text{max}}). It is this minimum wavelength, in conjuction with the spatial complexity of the model itself, which places the most restrictions on how we should build and discretize our model. No matter which numerical method we choose to solve the wave-equation we must have at least a certain number of "points" per minimum wavelength to ensure an accurate solution. In finite-difference (FD) simulations, these points are the discrete points of the FD grid. In spectral-element simulations (such as those performed by Salvus) these points are the GLL point locations with each spectral-element.
In Salvus, to properly balance accuracy and performance, we suggest to use anywhere from 6 - 9 points per minimum wavelength. When using standard 4th order elements (which have 5 GLL points along each edge), this equates to using 1.25 - 2 elements per wavelength. With this discussion in mind, one must first determine what frequency band they are interested in before one proceeds to the meshing stage, and this is why we now begin with the definition of our source wavelet. Below we choose a Ricker wavelet with a peak frequency of 5Hz, and plot the source-time function as well as the wavelet's frequency content.