Note

This page was generated from user_guide/dualLith/dualLithology.ipynb.
The binder version is quite slow due to memory and CPUs limitations
accordingly provided examples have been slightly simplified…
Interactive online version: Binder badge

Dual-lithology example

This example will run in approximately 20 minutes.

To run the notebook several post and pre-processing librairies are required:

  • panel

  • pyvista

  • pyevtk

  • jigsaw

To install these dependencies read the documentation in the user guide.

In addition, to these libraries several Python functions have been defined to ease the creation of gospl input files and the visualisation of the ouputs. These functions are located in the script folder in the same directory as this notebook:

  • buildMesh.py creates the initial mesh using the jigsaw library

  • runModel.py runs gospl based on the provided input file

  • readOutput.py builds a VTK file from model outputs

  • stratal.py extract the recorded stratigraphy in a specific region

[1]:
import numpy as np
import pyvista as pv
from script import stratal as strat
from script import readOutput as rout

import matplotlib
from matplotlib import cm
import matplotlib.pyplot as plt
import matplotlib.colors as colors
from mpl_toolkits.axes_grid1 import make_axes_locatable

label_size = 7
matplotlib.rcParams['xtick.labelsize'] = label_size
matplotlib.rcParams['ytick.labelsize'] = label_size
matplotlib.rc('font', size=6)

%matplotlib inline
%config InlineBackend.figure_format = 'svg'

1. Input files: topography and rainfall

First we have a data folder containing relevant inputs for the simulation.

  1. First, we will use a 0.1 degree resolution netcdf topography grid at present day (freely available paleoDEM files can be obtained from different sources, an example being the gPlates webportal), and

  2. Second, we will use a rainfall map based on current day precipitation, the dataset is again a netcdf file at 1.0 degree.

NOTE The folder is organised with specific sub-directory that are used later on during the mesh creation. You will need to follow this structure if you want to use your own dataset. Specifically the following folder names will be search: paleomap and precipitation.

2. Initial unstructured mesh generation

We start by creating the unstructured mesh using jigsaw library. It is ran from the terminal using buildMesh.py script.

python3 script/buildMesh.py -t=0 -d=data -s=5000,1000,100

This function takes 3 arguments:

  • t the time interval in Ma of the starting simulation time (here 0 Ma as we start at present day),

  • d the folder containing the input conditions. The code assumes that the paleo-elevations are located in the folder data under paleomap and are netCDF files of the form XXMa.nc with XX the specified time. It also assumes that the displacement maps are located in velocity and are off the form velocity_ XXMa.xy. Lastly for the paleo-precipitation, the assumed file is under precipitation and are netCDF files of the form XXMa.nc as well.

  • s is the space conditions for the jigsaw algorithm and consists of 3 values: the spacing in km for the mesh in the deep ocean (<=-1000 m), the spacing in km across shelf margin (>=-1000 m and < 0m) and the spacing in km in the continental domain. Here, the resolution will vary from 5000 km to 1000 km and 100 km respectively. This is set at pretty low resolution to ensure a quicker runtime. We often use the following values -s=100,30,15 in our run.

This function will take approximtely 5 minutes to complete…

[2]:
%run script/buildMesh.py -t=0 -d=data -s=5000,1000,100
Read scotese map (0.01 seconds)
Writing Jigsaw topo input .msh file (15.27 seconds)
Load topography grid (1.92 seconds)
Build space function (6.85 seconds)
Perform triangulation (22.71 seconds)
Get unstructured mesh (0.98 seconds)
Creation goSPL files took (48.55 seconds)

From the previous script, you will get a new folder created called input0 that contains the mesh information and the associated discretised elevation and rainfall values. Specifically you will have:

  • input0/mesh0.vtk: the VTK file containing the elevation at specified resolution,

  • input0/0Ma.npz: a Numpy compressed file with the mesh coordinates, cells, each vertice neighbours indices and the nodes elevation,

  • input0/rain0Ma.npz: a Numpy compressed file with the precipitation values for each mesh vertice.

Before going further, you can check the mesh validity by loading the created VTK file (input0/mesh0.vtk).

[3]:
mesh = pv.read('input0/mesh0.vtk')
elev = mesh.get_array(name='value')

earthRadius = 6.371e6
scale = 20.
factor = 1.+ (elev/earthRadius)*scale

mesh.points[:, 0] *= factor
mesh.points[:, 1] *= factor
mesh.points[:, 2] *= factor

contour = mesh.contour([0])

plot = pv.PlotterITK()
plot.add_mesh(mesh, scalars="value")
plot.add_mesh(contour, color="black", opacity=1.)
plot.show()

637a7f550c024f5fa334b9b6e8edaa80

3. Creating the initial stratigraphic mesh

In addition to the previous input files, we will specify an initial stratigraphic mesh. Here we impose a simple uniform 50 km thick layer for the entire globe composed of a 30% of fine and 70% of coarse sediment.

NOTE More complex stratigraphy can be defined with multiple layers varying spatially in thickness and composition.

In addition, to the layer thicknesses and compositions, gospl requires the porosity of each sediment type present in a given layer to be specified.

Here we follow the approach presented in the technical guide of the documentation, and the porosity is considered to varies with depth z.

\[\phi(z) = \phi_0 e^{−z/z_0}\]
[4]:
# Loading gospl mesh
loadMesh = np.load('input0/0Ma.npz')
gCoords = loadMesh["v"]
gZ = loadMesh["z"]

# Define layers variables
H = np.zeros((len(gZ),2)) # thickness
Z = H.copy()              # elevation
Fperc = H.copy()          # fine fraction
Fphi = H.copy()           # fine porosity
Sphi = H.copy()           # coarse porosity

H[:,1] = 50.0e3           # 50 km thick
Fperc[:,1] = 0.3          # 30% of fines
Z[:,0] = gZ - 50.e3       # elevation at the layer centre (base layer is assumed to have a 0 thickness)
Z[:,1] = gZ - 25.e3       # elevation at the layer centre (top layer centre elevation is at 25 km depth)

phis = 0.49               # Coarse sediment surface porosity
phif = 0.63               # Fine sediment surface porosity
z0s = 3700.0              # e-folding fine sediment thickness for porosity reduction 3700 m
z0f = 1960.0              # e-folding fine sediment thickness for porosity reduction 1960 m


# Compute porosity based on above equation
Fphi1 = phif * np.exp(-25.e3/z0f)
Sphi1 = phis * np.exp(-25.e3/z0s)
Fphi[:,1] = Fphi1
Sphi[:,1] = Sphi1

# Save the stratigraphic mesh as a Numpy file...
np.savez_compressed('input0/sed0Ma', strataH=H, strataF=Fperc, strataZ=Z, phiF=Fphi, phiS=Sphi)

4. Run gospl

Running gospl is done by calling the runModel.py script with the name of the input file as argument.

The Python script takes the following arguments:

  1. -i XXXX.yml specifying the input file name (required)

  2. -l True/False for outputing PETSc log during runtime (default is set to False)

  3. -v True/False for verbosing option during runtime (default is set to False)

You can open the input.yml file to look at the parameters that are setup for this model. A complete list of the gospl input variables is available in the user guide documentation.

This function will take approximately 5 minutes to complete with one CPU…

[5]:
# On a single processor...
%run script/runModel.py -i input.yml

# In parallel...
#!mpirun -np 4 python3 script/runModel.py -i input.yml
--- Initialisation Phase (0.85 seconds)
+++ Output Simulation Time: 0.00 years
--- Computational Step                       (2.28 seconds)
--- Computational Step                       (1.90 seconds)
--- Computational Step                       (2.15 seconds)
--- Computational Step                       (2.18 seconds)
--- Computational Step                       (1.93 seconds)
--- Computational Step                       (1.96 seconds)
--- Computational Step                       (2.29 seconds)
--- Computational Step                       (2.30 seconds)
--- Computational Step                       (2.18 seconds)
--- Computational Step                       (2.34 seconds)
+++ Output Simulation Time: 100000.00 years
--- Computational Step                       (2.50 seconds)
--- Computational Step                       (2.15 seconds)
--- Computational Step                       (2.20 seconds)
--- Computational Step                       (2.24 seconds)
--- Computational Step                       (2.16 seconds)
--- Computational Step                       (2.22 seconds)
--- Computational Step                       (2.26 seconds)
--- Computational Step                       (2.29 seconds)
--- Computational Step                       (2.44 seconds)
--- Computational Step                       (2.32 seconds)
+++ Output Simulation Time: 200000.00 years
--- Computational Step                       (2.71 seconds)
--- Computational Step                       (2.30 seconds)
--- Computational Step                       (2.32 seconds)
--- Computational Step                       (2.32 seconds)
--- Computational Step                       (2.32 seconds)
--- Computational Step                       (2.23 seconds)
--- Computational Step                       (2.32 seconds)
--- Computational Step                       (2.42 seconds)
--- Computational Step                       (2.34 seconds)
--- Computational Step                       (2.28 seconds)
+++ Output Simulation Time: 300000.00 years
--- Computational Step                       (2.89 seconds)
--- Computational Step                       (2.43 seconds)
--- Computational Step                       (2.51 seconds)
--- Computational Step                       (2.40 seconds)
--- Computational Step                       (2.46 seconds)
--- Computational Step                       (2.33 seconds)
--- Computational Step                       (2.50 seconds)
--- Computational Step                       (2.38 seconds)
--- Computational Step                       (2.46 seconds)
--- Computational Step                       (2.37 seconds)
+++ Output Simulation Time: 400000.00 years
--- Computational Step                       (3.14 seconds)
--- Computational Step                       (2.34 seconds)
--- Computational Step                       (2.60 seconds)
--- Computational Step                       (2.37 seconds)
--- Computational Step                       (2.55 seconds)
--- Computational Step                       (2.39 seconds)
--- Computational Step                       (2.47 seconds)
--- Computational Step                       (2.49 seconds)
--- Computational Step                       (2.50 seconds)
--- Computational Step                       (2.46 seconds)
+++ Output Simulation Time: 500000.00 years
--- Computational Step                       (3.30 seconds)
--- Computational Step                       (2.54 seconds)
--- Computational Step                       (2.56 seconds)
--- Computational Step                       (2.54 seconds)
--- Computational Step                       (2.56 seconds)
--- Computational Step                       (2.71 seconds)
--- Computational Step                       (2.71 seconds)
--- Computational Step                       (2.78 seconds)
--- Computational Step                       (2.81 seconds)
--- Computational Step                       (2.82 seconds)
+++ Output Simulation Time: 600000.00 years
--- Computational Step                       (3.73 seconds)
--- Computational Step                       (2.85 seconds)
--- Computational Step                       (2.83 seconds)
--- Computational Step                       (2.76 seconds)
--- Computational Step                       (2.80 seconds)
--- Computational Step                       (2.83 seconds)
--- Computational Step                       (2.94 seconds)
--- Computational Step                       (2.80 seconds)
--- Computational Step                       (2.99 seconds)
--- Computational Step                       (2.58 seconds)
+++ Output Simulation Time: 700000.00 years
--- Computational Step                       (4.07 seconds)
--- Computational Step                       (2.93 seconds)
--- Computational Step                       (3.06 seconds)
--- Computational Step                       (2.69 seconds)
--- Computational Step                       (3.09 seconds)
--- Computational Step                       (2.70 seconds)
--- Computational Step                       (3.11 seconds)
--- Computational Step                       (2.96 seconds)
--- Computational Step                       (3.02 seconds)
--- Computational Step                       (2.89 seconds)
+++ Output Simulation Time: 800000.00 years
--- Computational Step                       (4.13 seconds)
--- Computational Step                       (2.91 seconds)
--- Computational Step                       (3.04 seconds)
--- Computational Step                       (2.72 seconds)
--- Computational Step                       (2.94 seconds)
--- Computational Step                       (2.63 seconds)
--- Computational Step                       (2.94 seconds)
--- Computational Step                       (2.64 seconds)
--- Computational Step                       (2.96 seconds)
--- Computational Step                       (2.65 seconds)
+++ Output Simulation Time: 900000.00 years
--- Computational Step                       (4.18 seconds)
--- Computational Step                       (2.86 seconds)
--- Computational Step                       (2.94 seconds)
--- Computational Step                       (2.89 seconds)
--- Computational Step                       (3.04 seconds)
--- Computational Step                       (2.88 seconds)
--- Computational Step                       (2.96 seconds)
--- Computational Step                       (2.90 seconds)
--- Computational Step                       (3.21 seconds)
--- Computational Step                       (2.93 seconds)
+++ Output Simulation Time: 1000000.00 years
--- Computational Step                       (4.53 seconds)

+++
+++ Total run time (270.56 seconds)
+++

5. Visualisation in a notebook environment

The preferred way for visualising the model output is via Paraview by loading the time series file called gospl.xdmf available in the output folder (here called dual-lithologies).

Amongst the temporal variables outputed by gospl you will find:

  • surface elevation elev.

  • cumulative erosion & deposition values erodep.

  • flow accumulation flowAcc before pit filling.

  • flow accumulation fillAcc for depressionless surface.

  • river sediment load sedLoad.

  • fine sediment load sedLoadf when dual lithologies are accounted for.

  • uplift subsidence values if vertical tectonic forcing is considered uplift.

  • horizontal displacement values when considered hdisp.

  • precipitation maps based on forcing conditions rain.

Several filters, rendering and calculation can be done with Paraview but are beyond the scope of this example.

Here you will use the readOutput.py functions available in the script folder to visualise directly the model output in the notebook at the final time step.

The function requires several arguments:

  • path: the path to the input file

  • filename: the name of the input file

  • step: the step you wish to output (here set to 10 corresponding to the last output based on the input parameters: start time 0 year, end time 1 million years with an output every 0.1 million years)

  • back: set to False as the simulation is not a backward forward model

  • uplift: set to False as we are not considering any tectonic forcing

[6]:
# Reading the final output generated by gospl
output = rout.readOutput(path='./', filename='input.yml',step=10, back=False, uplift=False)

# Exporting the final output as a VTK mesh
output.exportVTK('step10.vtk')
Writing VTK file step10.vtk

We can now visualise the VTK output in the notebook directly:

[7]:
mesh = pv.read('step10.vtk')
elev = mesh.get_array(name='elev')

earthRadius = 6.371e6
scale = 20.
factor = 1.+ (elev/earthRadius)*scale

mesh.points[:, 0] *= factor
mesh.points[:, 1] *= factor
mesh.points[:, 2] *= factor

contour = mesh.contour([0])

plot = pv.PlotterITK()
plot.add_mesh(mesh, scalars="elev")
plot.add_mesh(contour, color="black", opacity=1.)
plot.show()

39f18e5f6dbf4e8ba3ca023a9ce46cc2

6. Extract stratigraphy

Finally, we will look at the recorded stratigraphy. The stratigraphic layer are recorded in gospl as HDF5 files stored in the output folder as stratal.XX.pX.h5 where XX is the output step and X the processor number.

The following information are stored:

  • elevation at time of deposition, considered to be to the current elevation for the top stratigraphic layer stratZ.

  • thickness of each stratigrapic layer stratH accounting for both erosion & deposition events.

  • proportion of fine sediment stratF contains in each stratigraphic layer.

  • porosity of coarse sediment phiS in each stratigraphic layer computed at center of each layer.

  • porosity of fine sediment phiF in each stratigraphic layer computed at center of each layer.

We will use the stratal.py function to extract the information above. It requires the following arguments: 1. path: the path to the input file 2. filename: the name of the input file 3. layer: the stratal file you wish to output

[8]:
# Load the function and specify the input file
strati = strat.stratal(path='./', filename='input.yml', layer=10)

# Read the stratigraphic dataset
strati.readStratalData()

# Interpolate the spherical dataset on a lon/lat regular grid
# by specifying the desired resolution and interpolation neighbours
strati.buildLonLatMesh(res=0.1, nghb=3)
Maximum layer number: 10
Current number of sedimentary: 52
Start building regular stratigraphic arrays
Percentage of arrays built : [####################] 100.0% DONE

We can visualise the create maps directly by doing…

[9]:
fig = plt.figure(figsize=(10,8))
ax = fig.add_subplot(111)

ax.imshow(np.flipud(strati.zi[-1,:,:]), extent=(-180, 180, -90, 90), vmin=-8000, vmax=8000, cmap=cm.bwr)
ax.set(xlabel='Longitude', ylabel='Latitude', yticks=np.arange(-90,120,30), xticks=np.arange(-180,180,30))
ax.minorticks_on()
../../_images/user_guide_dualLith_dualLithology_20_0.svg

We will now extract the stratigraphic layer for a specific region by using the writeMesh function which takes the following argument:

  • vtkfile the output name of the VTK stratigraphic mesh to create

  • lats latitude of the lower left and upper right corner of the region (specified between -90 and 90 degree)

  • lons longitude of the lower left and upper right corner of the region (specified between -180 and 180 degree)

  • sigma the standard deviation for Gaussian kernel as defined in SciPy gaussian_filter function.

The function returns the domain length in meters along the X and Y borders.

Here we chose a region around the Ganges–Brahmaputra–Meghna delta.

This function will take approximately 2 minutes to complete.

[10]:
length = strati.writeMesh(vtkfile='strati2D',
                          lats=[9,24],
                          lons=[82,97],
                          sigma=2.)

The function will build a VTK structured mesh containing the stratigraphic record for the region.

Here we will set the slice at the center of the domain…

We can viusualise the stratigraphic layers in the notebook:

[11]:
mesh = pv.read('strati2D.vts')
mesh.set_active_scalars('layID')
threshold = mesh.threshold([2,52])

surface = mesh.threshold([50,52])

# Position cross-section at the center of the region
slices = threshold.slice_orthogonal(x=length[0]/2, y=length[-1]/2, z=-10)

scale_factor = 5
slices[0].points[:, -1] *= scale_factor
slices[1].points[:, -1] *= scale_factor

contours0 = slices[0].contour(np.linspace(0, 51, 52))
contours1 = slices[1].contour(np.linspace(0, 51, 52))

surface.points[:, -1] *= scale_factor

p = pv.PlotterITK()
p.add_mesh(surface, scalars="dep elev", opacity=0.25)
p.add_mesh(slices[0], scalars="percfine")
p.add_mesh(slices[1], scalars="percfine")
p.add_mesh(contours0, color="black", opacity=1.)
p.add_mesh(contours1, color="black", opacity=1.)

p.show()

40332412656b49bdb7b59b14c4c913c0

[ ]: