```
# NBVAL_IGNORE_OUTPUT
"""
!apt-get update
!apt-get -qq install xvfb
!pip install pyvirtualdisplay
from pyvirtualdisplay import Display
display = Display(visible=0, size=(600, 400))
display.start()
""";
```

# 14 - Introduction to creating seismic synthetics

In this tutorial, we will create a simple geological model using GemPy, turn this into a 3D synthetic model of seismic wave velocities, then forward propagate a source through the model. We will then plot the synthetic shot gathers generated, and visualise the wavefield in 3D using PyVista.

## Side note (Google Colab, WLS, etc)

When running on some platforms (rule of thumb is those where creating floating windows is not supported), it is necessary to uncomment and run the following lines to prevent 3D renders from crashing the notebook. If you do not want to install these additional dependencies, then it is recommended that you comment out the GemPy and PyVista plotting cells

## Overview and setup

The synthetics which we will be building in this tutorial will be made with the use of GemPy, an open-source 3D geological modelling package for python. As this is not a core dependency of Devito, we will need to install it. If issues are encountered whilst installing GemPy into a `conda`

environment using `pip`

, you can alternatively create a python `venv`

and install Devito in this environment using `pip`

as per usual. Note that it will also be necesary to install an `ipykernel`

in this environment to run this notebook. From here, we can install GemPy:

```
# NBVAL_IGNORE_OUTPUT
try:
# Import gempy
import gempy as gp
except ModuleNotFoundError:
# Need to install these
! pip install aiohttp==3.7.1
! pip install pyvista==0.29
! pip install pyqt5
# Install gempy
! pip install gempy==2.2.9
# Import gempy
import gempy as gp
try:
# Import jinja2 (used for colour coding geology)
import jinja2
except ModuleNotFoundError:
# Install jinja2
! pip install jinja2
# Import jinja2
import jinja2
try:
# Check vtk notebook backend is installed
import ipyvtklink
except ModuleNotFoundError:
! pip install ipyvtklink
import ipyvtklink
```

`No module named 'osgeo'`

The simple geological model which we will be building is designed to evoke carbon-capture and storage (CCS) scenarios.

The model consists of a CO2 lens in a sandstone reservoir, with a shale layer in the overarching anticline providing the structural trap. This is then overlain by a layer of sediment, with water at the top of the model. Geological strata and their respective velocities are based on values detailed in Queißer et al. 2013, a paper imaging the P-wave velocity anomaly generated by CO2 injection into the Utsira Sand at Sleipnir in the North Sea using FWI. The model we will create features a similar shale trap/permeable sandstone reservoir structure, albeit with a small number of thick layers rather than the thin interbedding, to limit model complexity for this tutorial. Further inspiration was taken from Chadwick et al. 2004, a paper characterizing the Utsira Sand reservoir based on 2D seismic lines and well logs.

## Creating our geological model:

To begin, alongside GemPy, we need to import auxiliary modules:

```
# Import auxiliary modules
import numpy as np
%matplotlib inline
```

We will now set up a GemPy `Model`

object. This encapsulates the grid onto which the scalar fields associated with various surfaces are interpolated. Note that the extent is slightly greater than it will be for our Devito model (an extra half a grid spacing is added to each side).

A comparison of the cell-centered vs node-centered conventions of GemPy and Devito respectively, along with the differences in how they measure extent. It is necessary to account for this to ensure that the two grids are co-located.

As we can see in the figure above, this is due to differences in the way in which grids are defined in each package and is necessary to ensure that the model is not stretched and distorted when transistioning between the two, and that they are correctly aligned.

```
# Set overarching model parameters
= (-5., 1005., -5., 1005., -1005., 5.)
extent = (101, 101, 101)
shape
= gp.create_model('Gempy-tutorial')
geo_model = gp.init_data(geo_model, extent=extent, resolution=shape) geo_model
```

`Active grids: ['regular']`

Setting up Theano for our model (used by GemPy for interpolation). Bear in mind that that this may take some time to run.

```
# NBVAL_IGNORE_OUTPUT
=['geology'], theano_optimizer='fast_compile') gp.set_interpolator(geo_model, output
```

```
Setting kriging parameters to their default values.
Compiling theano function...
Level of Optimization: fast_compile
Device: cpu
Precision: float64
Number of faults: 0
Compilation Done!
Kriging values:
values
range 1749.371316
$C_o$ 72864.285714
drift equations [3]
```

`<gempy.core.interpolator.InterpolatorModel at 0x7f1ddc20fe50>`

As the top CO2 surface is truncated by the upper shale layer, we will need to separate the geological strata into two GemPy `Series`

. Each `Series`

object, as the name suggests is intended to correspond with a geological unit, and they can be made to onlap, erode, etc one another. Whilst in practice, the top CO2 contact is not an erosive surface, treating it as such is the most straightforward way to create the desired truncation.

A default series is included in the model. As such, rather than creating a new series, we will simply rename it to ‘Lower’. As you can imagine, this is going to be used to contain the lower geological units, these being the lower shale, reservoir sandstone, and CO2 lens.

`'Default series': 'Lower'}) geo_model.rename_series({`

And now add our surfaces:

```
# NBVAL_IGNORE_OUTPUT
'co2', 'sands', 'lowershale']) geo_model.add_surfaces([
```

surface | series | order_surfaces | color | id | |
---|---|---|---|---|---|

0 | co2 | Lower | 1 | #015482 | 1 |

1 | sands | Lower | 2 | #9f0052 | 2 |

2 | lowershale | Lower | 3 | #ffbe00 | 3 |

We will now set some points for the base of the sands and CO2. The lower shale is considered the basement, meaning that its base does not need to be defined and it will extend to the bottom of the model. Alongside these points, we wil need to define an orientation for the surface.

To minimise repetition, we will define a function to loop over a list of points and add each to the surface.

```
def create_surface(model, points, surface):
"""Add a list of points to a surface in a model"""
= ('X', 'Y', 'Z')
xyz for point in points:
= {**dict(zip(xyz, point)), 'surface': surface}
kwargs **kwargs)
model.add_surface_points(
# The points defining the base of the sand layer
= [(322, 135, -783), (635, 702, -791), (221, 668, -772), (732, 235, -801), (442, 454, -702)]
sand_points
# Call our function
'sands')
create_surface(geo_model, sand_points,
# Add the surface orientation
=442., Y=495., Z=-752.,
geo_model.add_orientations(X='sands', pole_vector=(0.05, 0.05, 0.95)) surface
```

X | Y | Z | G_x | G_y | G_z | smooth | surface | |
---|---|---|---|---|---|---|---|---|

0 | 442.0 | 495.0 | -752.0 | 0.05 | 0.05 | 0.95 | 0.01 | sands |

We will now repeat this process for the CO2 lens.

```
# Points defining the base of the CO2 layer
= [(322, 135, -650), (635, 702, -650), (221, 668, -650), (732, 235, -650), (442, 454, -650)]
co2_points
'co2')
create_surface(geo_model, co2_points,
# Add the surface orientation
=495., Y=495., Z=-650.,
geo_model.add_orientations(X='co2', pole_vector=(0., 0., 1.)) surface
```

X | Y | Z | G_x | G_y | G_z | smooth | surface | |
---|---|---|---|---|---|---|---|---|

1 | 495.0 | 495.0 | -650.0 | 0.00 | 0.00 | 1.00 | 0.01 | co2 |

0 | 442.0 | 495.0 | -752.0 | 0.05 | 0.05 | 0.95 | 0.01 | sands |

We will now add an upper series, containing statigraphy above the CO2 lens.

`'Upper') geo_model.add_series(`

order_series | BottomRelation | isActive | isFault | isFinite | |
---|---|---|---|---|---|

Lower | 1 | Erosion | True | False | False |

Upper | 2 | Erosion | False | False | False |

As we can see, the upper series has been added below the lower series. This is not ideal for obvious reasons, and hence we will reorder them:

`'Upper', 'Lower']) geo_model.reorder_series([`

order_series | BottomRelation | isActive | isFault | isFinite | |
---|---|---|---|---|---|

Upper | 1 | Erosion | False | False | False |

Lower | 2 | Erosion | True | False | False |

And add our remaining surfaces:

```
# NBVAL_IGNORE_OUTPUT
'water', 'sediments', 'uppershale']) geo_model.add_surfaces([
```

surface | series | order_surfaces | color | id | |
---|---|---|---|---|---|

0 | co2 | Lower | 1 | #015482 | 1 |

1 | sands | Lower | 2 | #9f0052 | 2 |

2 | lowershale | Lower | 3 | #ffbe00 | 3 |

3 | water | Lower | 4 | #728f02 | 4 |

4 | sediments | Lower | 5 | #443988 | 5 |

5 | uppershale | Lower | 6 | #ff3f20 | 6 |

As these surfaces are not mapped to the upper series by default, we shall do so:

```
# NBVAL_IGNORE_OUTPUT
'Upper': ('water', 'sediments', 'uppershale')}) gp.map_stack_to_surfaces(geo_model, {
```

surface | series | order_surfaces | color | id | |
---|---|---|---|---|---|

3 | water | Upper | 1 | #728f02 | 1 |

4 | sediments | Upper | 2 | #443988 | 2 |

5 | uppershale | Upper | 3 | #ff3f20 | 3 |

0 | co2 | Lower | 1 | #015482 | 4 |

1 | sands | Lower | 2 | #9f0052 | 5 |

2 | lowershale | Lower | 3 | #ffbe00 | 6 |

Now we will add the points for the upper series. Note that there is only a single orientation included. It is not necessary to define an orientation for each surface, so long as there is an orientation in the series.

```
# Surface points
= [(322, 135, -633), (635, 702, -641), (221, 668, -622), (732, 235, -651), (442, 454, -552)]
uppershale_points = [(322, 135, -433), (635, 702, -441), (221, 668, -422), (732, 235, -451), (442, 454, -352)]
sediments_points = [(232, 153, -221), (653, 234, -216), (112, 872, -198), (532, 572, -223),
water_points 722, 884, -189), (632, 429, -201), (732, 348, -222)]
(
# Add the points to our surfaces
'uppershale')
create_surface(geo_model, uppershale_points, 'sediments')
create_surface(geo_model, sediments_points, 'water')
create_surface(geo_model, water_points,
# Set an orientation
=442., Y=495., Z=-502.,
geo_model.add_orientations(X='uppershale', pole_vector=(0.05, 0.05, 0.95)) surface
```

X | Y | Z | G_x | G_y | G_z | smooth | surface | |
---|---|---|---|---|---|---|---|---|

2 | 442.0 | 495.0 | -502.0 | 0.05 | 0.05 | 0.95 | 0.01 | uppershale |

1 | 495.0 | 495.0 | -650.0 | 0.00 | 0.00 | 1.00 | 0.01 | co2 |

0 | 442.0 | 495.0 | -752.0 | 0.05 | 0.05 | 0.95 | 0.01 | sands |

Finally, we can add the p wave velocities associated with each of these layers. Note that any parameter can be set in this manner (density, elastic parameters, attenuation, etc) if desired for more complex synthetics.

```
# NBVAL_IGNORE_OUTPUT
1.5, 1.75, 2.5, 1.1, 2., 2.5], ['vp'])
geo_model.add_surface_values([ geo_model.surfaces
```

surface | series | order_surfaces | color | id | vp | |
---|---|---|---|---|---|---|

3 | water | Upper | 1 | #728f02 | 1 | 1.500000 |

4 | sediments | Upper | 2 | #443988 | 2 | 1.750000 |

5 | uppershale | Upper | 3 | #ff3f20 | 3 | 2.500000 |

0 | co2 | Lower | 1 | #015482 | 4 | 1.100000 |

1 | sands | Lower | 2 | #9f0052 | 5 | 2.000000 |

2 | lowershale | Lower | 3 | #ffbe00 | 6 | 2.500000 |

Now we can visualise our model, plotting data points and orientations. Then we must compute our model to interpolate the surfaces and any associated scalar fields. Then we can plot our surfaces and the associated units.

```
# NBVAL_IGNORE_OUTPUT
# Compute the model. Note that a solution is returned. We will use this later
= gp.compute_model(geo_model) sol
```

```
# NBVAL_SKIP
# Set up plotter
= gp.plot_3d(geo_model, plotter_type='background', notebook=True)
p3d # Plot data points and orientations
p3d.plot_data()
# Plot the surfaces
p3d.plot_surfaces()# Plot the lithological units
'lith') p3d.plot_structured_grid(
```

```
[StructuredGrid (0x7f1d171d7360)
N Cells: 1000000
N Points: 1030301
X Bounds: 0.000e+00, 1.000e+03
Y Bounds: 0.000e+00, 1.000e+03
Z Bounds: -1.000e+03, 0.000e+00
Dimensions: 101, 101, 101
N Arrays: 1,
<matplotlib.colors.ListedColormap at 0x7f1d5bb3a8d0>]
```

## Bridging the gap from GemPy to Devito:

As you may have noticed, when we compute our GemPy model, a `Solution`

object is returned. From this, we can extract the rasterized values attached to each of our geological units. With this in mind, we can print the solution values:

` sol.values_matrix`

```
array([[2.5 , 2.44720385, 2. , ..., 1.5 , 1.5 ,
1.5 ]])
```

You will notice that these values correspond with the p wave velocities we specified. However, they are in the form of 1D vector. Consequently, will need to reshape this array to fit into the `vp`

parameter of a Devito `Model`

. Note that you could do this with further parameters such as density or shear wave velocity for more complex models. In this case, you would want to set up a Devito `Function`

to contain each parameter.

Note that in this case, we need to select c-like index order to get the axis in the correct order.

```
# Reshaping our data to the shape required by Devito
= np.reshape(sol.values_matrix, shape, order='C')
reshaped reshaped.shape
```

`(101, 101, 101)`

Now let us plot a slice through this model for quality checking purposes.

```
# NBVAL_IGNORE_OUTPUT
import matplotlib.pyplot as plt
# Take the center slice in the x direction
# Remember that in Devito, indexing convention is [x, y, z] (need to flip for correct imshow display)
50].T, cmap='viridis', origin='lower')
plt.imshow(reshaped[
plt.colorbar() plt.show()
```

## Seismic modelling with Devito

We can now start building our Devito model. The following draws heavily from the Devito `examples/sesimic/tutorials/01_modelling.ipynb`

notebook. We will begin, as always with some imports.

```
import devito as dv
from examples.seismic import Model
```

As mentioned earlier, Devito and GemPy have slightly different grid implementations, so we need to tweak the Devito configuration slightly to make it map to the GemPy grid. We can now construct a Devito `Model`

. This is a convenience object encapsulating the necessary parameters and components of an acoustic wave model, including additional damping layers around the perimeter (specified by `bcs="damp"`

). For custom setups, see the `examples/userapi/04_boundary_conditions.ipynb`

. Note that we are using a relatively large amount of damping layers here. This is to avoid our gathers becoming too messy, and ensure that reflections from horizons can be straightforwardly identified in the gathers.

```
# NBVAL_IGNORE_OUTPUT
= Model(vp=reshaped, origin=(0., 0., -1000.), spacing=(10., 10., 10.), shape=shape, nbl=30, space_order=4, bcs="damp") seis_model
```

```
Operator `initdamp` ran in 0.02 s
Operator `pad_vp` ran in 0.02 s
```

Now we will set up the time axis for our model. Again, this is a convenience object, which we will use in setting up the source and recievers.

```
from examples.seismic import TimeAxis
= 0. # Simulation starts a t=0
t0 = 1000. # Simulation last 1 second (1000 ms)
tn = seis_model.critical_dt # Time step from model grid spacing
dt
= TimeAxis(start=t0, stop=tn, step=dt) time_range
```

We will position our source at a depth of 20m, and center it in all other axes.

```
# NBVAL_IGNORE_OUTPUT
from examples.seismic import RickerSource
= 0.015 # Source peak frequency is 15Hz (0.015 kHz)
f0 = RickerSource(name='src', grid=seis_model.grid, f0=f0,
src =1, time_range=time_range)
npoint
# First, position source centrally in all dimensions, then set depth
= np.array(seis_model.domain_size) * .5
src.coordinates.data[:] 0, -1] = -20 # Depth is 20m
src.coordinates.data[
# We can plot the time signature to see the wavelet
src.show()
```

We will also configure our recievers in a line along the x axis, centered in the y, also at a depth of 20m.

```
from examples.seismic import Receiver
# Create symbol for 101 receivers
= Receiver(name='rec', grid=seis_model.grid, npoint=101, time_range=time_range)
rec
# Prescribe even spacing for receivers along the x-axis
0] = np.linspace(0, seis_model.domain_size[0], num=101)
rec.coordinates.data[:, 1] = 0.5*seis_model.domain_size[1]
rec.coordinates.data[:, -1] = -20. # Depth is 20m rec.coordinates.data[:,
```

In Devito, equation parameters which vary in space only are represented using `Function`

objects. If we also want them to vary over time, we must use a `TimeFunction`

.

With this, we can define our partial differential equation.

```
# Define the wavefield with the size of the model and the time dimension
= dv.TimeFunction(name="u", grid=seis_model.grid, time_order=2, space_order=4)
u
# We can now write the PDE
= seis_model.m * u.dt2 - u.laplace + seis_model.damp * u.dt
pde
# The PDE representation is as on paper
pde
```

\(\displaystyle \operatorname{damp}{\left(x,y,z \right)} \frac{\partial}{\partial t} u{\left(t,x,y,z \right)} - \frac{\partial^{2}}{\partial x^{2}} u{\left(t,x,y,z \right)} - \frac{\partial^{2}}{\partial y^{2}} u{\left(t,x,y,z \right)} - \frac{\partial^{2}}{\partial z^{2}} u{\left(t,x,y,z \right)} + \frac{\frac{\partial^{2}}{\partial t^{2}} u{\left(t,x,y,z \right)}}{\operatorname{vp}^{2}{\left(x,y,z \right)}}\)

Now create our update stencil:

```
# NBVAL_IGNORE_OUTPUT
# This discrete PDE can be solved in a time-marching way updating u(t+dt) from the previous time step
# Devito as a shortcut for u(t+dt) which is u.forward. We can then rewrite the PDE as
# a time marching updating equation known as a stencil using customized SymPy functions
= dv.Eq(u.forward, dv.solve(pde, u.forward))
stencil stencil
```

\(\displaystyle u{\left(t + dt,x,y,z \right)} = \frac{- \frac{- \frac{2.0 u{\left(t,x,y,z \right)}}{dt^{2}} + \frac{u{\left(t - dt,x,y,z \right)}}{dt^{2}}}{\operatorname{vp}^{2}{\left(x,y,z \right)}} + \frac{\partial^{2}}{\partial x^{2}} u{\left(t,x,y,z \right)} + \frac{\partial^{2}}{\partial y^{2}} u{\left(t,x,y,z \right)} + \frac{\partial^{2}}{\partial z^{2}} u{\left(t,x,y,z \right)} + \frac{\operatorname{damp}{\left(x,y,z \right)} u{\left(t,x,y,z \right)}}{dt}}{\frac{\operatorname{damp}{\left(x,y,z \right)}}{dt} + \frac{1}{dt^{2} \operatorname{vp}^{2}{\left(x,y,z \right)}}}\)

Now we can set up our source and reciever terms to include in our `Operator`

.

```
# Finally we define the source injection and receiver read function to generate the corresponding code
= src.inject(field=u.forward, expr=src * dt**2 / seis_model.m)
src_term
# Create interpolation expression for receivers
= rec.interpolate(expr=u.forward) rec_term
```

Create our operator:

`= dv.Operator([stencil] + src_term + rec_term, subs=seis_model.spacing_map) op `

And run it.

```
# NBVAL_IGNORE_OUTPUT
=time_range.num-1, dt=seis_model.critical_dt) op(time
```

`Operator `Kernel` ran in 8.26 s`

```
PerformanceSummary([(PerfKey(name='section0', rank=None),
PerfEntry(time=8.240707, gflopss=0.0, gpointss=0.0, oi=0.0, ops=0, itershapes=[])),
(PerfKey(name='section1', rank=None),
PerfEntry(time=0.0006800000000000077, gflopss=0.0, gpointss=0.0, oi=0.0, ops=0, itershapes=[])),
(PerfKey(name='section2', rank=None),
PerfEntry(time=0.011009999999999985, gflopss=0.0, gpointss=0.0, oi=0.0, ops=0, itershapes=[]))])
```

We can now plot our shot record using everyone’s favourite colourmap. We can clearly see the reflected arrivals from the seabed, top shale, and top CO2. We can also distinguish the base of the CO2 lens and the interface between the reservoir sandstone and the underlying shale.

```
# NBVAL_IGNORE_OUTPUT
='viridis', aspect='auto', vmax=0.01, vmin=-0.01)
plt.imshow(rec.data, cmap"Reciever number")
plt.xlabel("Time (ms)")
plt.ylabel(
plt.colorbar() plt.show()
```

## Visualisation with PyVista:

As PyVista is a dependency of GemPy, we can use its plotting capabilities to plot some slices to visualise our wavefield

```
import pyvista as pv
# Set default pyvista backend
'ipyvtklink')
pv.set_jupyter_backend(
# Trim down the data from u to remove damping field
= u.data[1, 30:-30, 30:-30, 30:-30]
trimmed_data
# Create the spatial reference
= pv.UniformGrid()
grid
# Set the grid dimensions: shape + 1 because we want to inject our values on
# the CELL data
= np.array(trimmed_data.shape) + 1
grid.dimensions
# Edit the spatial reference
= (0., 0., -1000.) # The bottom left corner of the data set
grid.origin = (10, 10, 10) # These are the cell sizes along each axis grid.spacing
```

We can now fill the grid cells:

```
# Add the data values to the cell data
"values"] = trimmed_data.flatten(order="F") # Flatten the array! grid.cell_data[
```

And plot!

```
# NBVAL_SKIP
= grid.slice_orthogonal(x=200, y=200, z=-500)
orth_slices
='seismic', clim=[-0.01, 0.01]) orth_slices.plot(cmap
```

```
[(2456.1701691039184, 2456.1701691039184, 1456.1701691039184),
(505.0, 505.0, -495.0),
(0.0, 0.0, 1.0)]
```

Some other visualisation shenanigans:

```
# NBVAL_SKIP
= grid.slice_along_axis(n=5, axis="y")
y_slices = pv.Plotter()
p ="k")
p.add_mesh(grid.outline(), color='seismic', clim=[-0.01, 0.01], opacity=0.8)
p.add_mesh(y_slices, cmap p.show()
```

```
[(2456.1701691039184, 2456.1701691039184, 1456.1701691039184),
(505.0, 505.0, -495.0),
(0.0, 0.0, 1.0)]
```