Inputs definition¶
For this second recipe, we will use the same inputs as the ones defined for Recipe 1.
Inputs definition
Used inputs are:
an initial mesh with the Gulf of Mexico as the region of interest
an elevation set to present-day toppography based on ETOPO5 dataset
a rainfall grid extracted from CPC collection of precipitation data sets containing global monthly values since 1979
a generic tectonic conditions (uplift/subsidence) where regions above 500 m are uplifted (0.1 cm/yr) and the ones below -10 m are subsiding (-0.1 cm/yr)
Tip
For building these forcings, you might want to go to the corresponding Jupyter notebook here.
Visualising inputs from Recipe 1¶
import meshio
import numpy as np
import pyvista as pv
paleovtk = 'gospl_data/mesh25km.vtk'
mesh = pv.read(paleovtk)
elev = mesh.get_array(name='elev')
earthRadius = 6.371e6
scale = 10.
factor = 1.+ (elev/earthRadius)*scale
mesh.points[:, 0] *= factor
mesh.points[:, 1] *= factor
mesh.points[:, 2] *= factor
contour = mesh.contour([0])
plotter = pv.PlotterITK()
plotter.add_mesh(mesh, scalars="elev")
plotter.add_mesh(contour, color="black", opacity=1.)
plotter.show()
Building sedimentary layers¶
Multi-lithologies stratigraphy¶
We will first create an inital stratigraphic architecture consisting of 3 layers.
It is possible to build a complex layering with variable thicknesses and proportions for each layer. Here, for the sake of simplicity, we will only define layers which have constant spatial characteristics.
Stratigraphic layers in gospl are defined on the unstructured mesh, so first we load this mesh…
# Loading gospl mesh
loadMesh = np.load('gospl_data/mesh25km.npz')
gCoords = loadMesh["v"]
gZ = loadMesh["z"]
Each layer corresponds to a specific time interval and is characterised by a number of parameters:
layer thickness (meters)
percentage of fine sediment
percentatage of weathered sediment
elevation at the time of deposition
averaged porosity for coarse sediment
averaged porosity for fine sediment
averaged porosity for weathered sediment
Let’s define each of these variables:
# Define layers variables
H = np.zeros((len(gZ),3)) # thickness
Z = H.copy() # elevation
Fperc = H.copy() # fine fraction
Wperc = H.copy() # weathered fraction
Fphi = H.copy() # fine porosity
Sphi = H.copy() # coarse porosity
Wphi = H.copy() # weathered porosity
Now we will set the initial thicknesses and compositions of the layers assuming uniform parameters:
H[:,0] = 1.0 # 1 m thick
H[:,1] = 5.0e3 # 10 km thick
H[:,2] = 10.0e3 # 10 km thick
Fperc[:,0] = 0.4 # 40% of fines
Fperc[:,1] = 0.4 # 40% of fines
Fperc[:,2] = 0.4 # 40% of fines
Wperc[:,0] = 0.2 # 20% of fines
Wperc[:,1] = 0.2 # 20% of fines
Wperc[:,2] = 0.2 # 20% of fines
Z[:,0] = gZ - 15000.5 # elevation at the centre of layer 0
Z[:,1] = gZ - 12500.0 # elevation at the centre of layer 1
Z[:,2] = gZ - 5000.0 # elevation at the centre of layer 2
To estimate the porosities we will assume that compaction is dependent on deposition rate and sediment surface porosity following the formulation proposed gospl and explains in more details here.
As such, we first define:
surface porosity of sediments,
the e-folding sediment thickness for porosity reduction.
phis = 0.49 # Coarse sediment surface porosity
phif = 0.63 # Fine sediment surface porosity
phiw = 0.65 # Weathered sediment surface porosity
z0s = 3700.0 # e-folding coarse sediment thickness for porosity reduction 3700 m
z0f = 1960.0 # e-folding fine sediment thickness for porosity reduction 1960 m
z0w = 1600.0 # e-folding weathered sediment thickness for porosity reduction 1960 m
Using Sclater and Christie, 1980 based on many sedimentary basins observations we then get:
# Weathered sediment porosity for each layer
Wphi0 = phiw * np.exp(-15000.5/z0w)
Wphi1 = phiw * np.exp(-12500./z0w)
Wphi2 = phiw * np.exp(-5000./z0w)
Wphi[:,0] = Wphi0
Wphi[:,1] = Wphi1
Wphi[:,2] = Wphi2
# Fine sediment porosity for each layer
Fphi0 = phif * np.exp(-15000.5/z0f)
Fphi1 = phif * np.exp(-12500./z0f)
Fphi2 = phif * np.exp(-5000./z0f)
Fphi[:,0] = Fphi0
Fphi[:,1] = Fphi1
Fphi[:,2] = Fphi2
# Coarse sediment porosity for each layer
Sphi0 = phis * np.exp(-15000.5/z0s)
Sphi1 = phis * np.exp(-12500./z0s)
Sphi2 = phis * np.exp(-5000./z0s)
Sphi[:,0] = Sphi0
Sphi[:,1] = Sphi1
Sphi[:,2] = Sphi2
We now save the stratigraphic grid (with all defined gospl parameters):
# Save the stratigraphic mesh as a Numpy file...
np.savez_compressed('gospl_data/sedlay25km', strataH=H, strataF=Fperc, strataW=Wperc, strataZ=Z,
phiF=Fphi, phiS=Sphi, phiW=Wphi)
Spatially variable layer composition¶
Note
To build a spatially variable layer composition one could for example extract specific subsurface shapefiles and define the layers accordingly.
Here we will define a surface layer with one side of the Mississippi drainage basin corresponding to one lithology and the other to a second one.
def xyz2lonlat(coords, radius=6378137.0):
"""
Convert x,y,z representation of cartesian points of the
spherical triangulation to lat/lon.
"""
gLonLat = np.zeros((len(coords), 2))
gLonLat[:, 1] = np.arcsin(coords[:, 2] / radius)
gLonLat[:, 0] = np.arctan2(coords[:, 1], coords[:, 0])
gLonLat[:, 1] = np.mod(np.degrees(gLonLat[:, 1]) + 90, 180.0) -90.
gLonLat[:, 0] = np.mod(np.degrees(gLonLat[:, 0]) + 180.0, 360.0) -180.
return gLonLat
# Loading gospl mesh
loadMesh = np.load('gospl_data/mesh25km.npz')
gCoords = loadMesh["v"]
gpoints = len(gCoords)
gZ = loadMesh["z"]
gCells = loadMesh["c"]
# Convert to lon/lat coordinates
gLonLat = xyz2lonlat(gCoords)
We now get the points East of longitude -95:
idSed1 = np.where(gLonLat[:,0]>-95)[0]
As for the previous example, we now need to specify the characteristics of each layer.
# Define layers variables
H = np.zeros((len(gZ),3)) # thickness
Z = H.copy() # elevation
Fperc = H.copy() # fine fraction
Wperc = H.copy() # weathered fraction
Fphi = H.copy() # fine porosity
Sphi = H.copy() # coarse porosity
Wphi = H.copy() # weathered porosity
H[:,0] = 1.0 # 1 m thick
H[:,1] = 1.0 # 1 m thick
H[:,2] = 10.0e3 # 10 km thick
We now define the proportion of each sediment present in the sedimentary layers:
# We set the points East of longitude -95 with a specific lithology
Fperc[idSed1,0] = 1.0 # 100% of fines
Fperc[idSed1,1] = 1.0 # 100% of fines
Fperc[idSed1,2] = 1.0 # 100% of fines
# We do not consider any weathered sediments in the stratigraphic layers
Wperc[:,0] = 0. # 0% of weathered
Wperc[:,1] = 0. # 0% of weathered
Wperc[:,2] = 0. # 0% of weathered
Here again, we estimate the porosities assuming that compaction is dependent on deposition rate and sediment surface porosity following the formulation proposed gospl:
Z[:,0] = gZ - 10001.5 # elevation at the centre of layer 0
Z[:,1] = gZ - 10000.5 # elevation at the centre of layer 1
Z[:,2] = gZ - 5000.0 # elevation at the centre of layer 2
phis = 0.49 # Coarse sediment surface porosity
phif = 0.49 # Fine sediment surface porosity
phiw = 0.65 # Weathered sediment surface porosity
z0s = 3700.0 # e-folding coarse sediment thickness for porosity reduction 3700 m
z0f = 3700.0 # e-folding fine sediment thickness for porosity reduction 1960 m
z0w = 1600.0 # e-folding weathered sediment thickness for porosity reduction 1960 m
# Compute porosity based on above equation
Wphi0 = phiw * np.exp(-10001.5/z0w)
Wphi1 = phiw * np.exp(-10000.5/z0w)
Wphi2 = phiw * np.exp(-5000./z0w)
Wphi[:,0] = Wphi0
Wphi[:,1] = Wphi1
Wphi[:,2] = Wphi2
Fphi0 = phif * np.exp(-10001.5/z0f)
Fphi1 = phif * np.exp(-10000.5/z0f)
Fphi2 = phif * np.exp(-5000./z0f)
Fphi[:,0] = Fphi0
Fphi[:,1] = Fphi1
Fphi[:,2] = Fphi2
Sphi0 = phis * np.exp(-10001.5/z0s)
Sphi1 = phis * np.exp(-10000.5/z0s)
Sphi2 = phis * np.exp(-5000./z0s)
Sphi[:,0] = Sphi0
Sphi[:,1] = Sphi1
Sphi[:,2] = Sphi2
We now save the stratigraphic grid
# Save the stratigraphic mesh as a Numpy file...
np.savez_compressed('gospl_data/surflay25km', strataH=H, strataF=Fperc, strataW=Wperc, strataZ=Z,
phiF=Fphi, phiS=Sphi, phiW=Wphi)
Let’s have a look at the surface sedimentary layer:
vtk_file = 'gospl_data/surflay25km.vtk'
gCoords = loadMesh["v"]
gpoints = len(gCoords)
gZ = loadMesh["z"]
gCells = loadMesh["c"]
vis_mesh = meshio.Mesh(gCoords, {"triangle": gCells},
point_data={"elev": gZ,
"propSed":Fperc[:,2] ,
})
meshio.write(vtk_file, vis_mesh)
mesh = pv.read(vtk_file)
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])
plotter = pv.PlotterITK()
plotter.add_mesh(mesh, scalars="propSed")
plotter.add_mesh(contour, color="black", opacity=1.)
plotter.show()