Analyse pyBadlands stratigraphic output

If the stratigraphic structure is turned on in the XmL input file, pyBadlands produces sedimentary layers Hdf5 files. The stratigraphic layers are defined on a regularly spaced grid and a layer is recorded at each layer time interval given by the user.

Here we show how we can visualise quickly the structure of the stratigraphic layer in an IPython notebook.



In [ ]:

    
import numpy as np
import matplotlib.mlab as ml
import matplotlib.pyplot as plt
import scipy.ndimage.filters as filters
from scipy.interpolate import RectBivariateSpline
from scipy.ndimage.filters import gaussian_filter

import pybadlands_companion.stratalAnalyse_DT as strata

# display plots in SVG format
%config InlineBackend.figure_format = 'svg'
%matplotlib inline

Loading stratigraphic file

First we need to load the stratigraphic file. The files are located in the h5/ folder in the simulation main output folder and are named using the following convention:

sed.timeT.pX.hdf5

with T the display time index and X the number of the partition (used in the parallel version). In cases where you ran your simulation in parallel you will also need to give the number of CPUs used (cpus).

To load a file you will need to give the folder name and the number of processors used in your simulation.

For more information regarding the function uncomment the following line.



In [ ]:

    
# help(strata.stratalSection.__init__)



In [ ]:

    
folder = '/workspace/volume/Examples/dyntopo1/output/h5/'
strat = strata.stratalSection(folder)

Then we need to load a particular output time interval (this is the T parameter in the hdf5 file name convention).

Note

This number is not always the number of sedimentary layers for this particular time step as you could have chosen in the input file to have more than 1 sedimentary layer recorded by output interval!



In [ ]:

    
# help(strat.loadStratigraphy)



In [ ]:

    
strat.loadStratigraphy(100)

Building a cross-section

We then slice the stratigraphic mesh to visualise the sedimentary architecture along a given cross-section.

To create the cross-section you will need to provide:

the position of the segment in the simulation space (xo,yo) and (xm,ym),
a gaussian filter value for smoothing (gfilt a value of 0 can be used for non-smoothing) and
the resolution of the cross-section (based on a number of points: nbpts).



In [ ]:

    
# Build a cross-section along X axis
x1 = 0.
x2 = 300000.
y1 = 100000
y2 = 100000

# Interpolation parameters
nbpts = 301
gfilt = 1



In [ ]:

    
# help(strat.buildSection)



In [ ]:

    
strat.buildSection(xo = x1, yo = y1, xm = x2, ym = y2, pts = nbpts, gfilter = gfilt)

Stratal stacking pattern

We use plotly to visualise the vertival cross-section of stratal stacking apttern.



In [ ]:

    
#help(strata.viewSection)



In [ ]:

    
strata.viewSection(width = 1000, height = 500, cs = strat, 
            dnlay = 2, rangeX=[220000, 300000], rangeY=[-1050,-50],
            linesize = 0.5, title='Stratal stacking pattern coloured by time')

Sequence stratigraphy methods proposed in this notebook

New analytical methods are proposed in recent years on the interpretation of depositional systems. Here, we apply three methods, including

(i) the systems tracts method in which the declaration of stratigraphic sequences is based on relative sea-level change;
(ii) the shoreline trajectory analysis that defines different trajectory classes according to the migration of shoreline;
(iii) the accommodation sucession method that focuses on the objective observation of depositional trends and then defines different sequence sets reponding to the competing ration between the change of accommodation and sediment supply.

1- Systems Tracts Method - based on relative sea-level

There are several models of systems tracts within depositional sequences, here we use the most simple one called the four systems tract model (figure is from Helland-Hansen & Hampson (2009)):

highstand systems tract HST
falling-stage systems tract FSST
lowstand systems tract LST
transgressive systems tract TST

For each of these system tracts we define a given color. We use the 'colorlover' library link

Relative base-level

We first visualize the sea-level curve, and then extract manually the time intervals that bound different systems tracts, as showed in the above figure.



In [ ]:

    
time, Sealevel = strata.readSea('data/sealevel.csv')
timeStep = 1e5

# Plot sea-level
strata.viewData(x0 = time/timeStep, y0 = Sealevel, width = 600, height = 400, linesize = 3, 
                color = '#6666FF',title='Sea-level curve',xlegend='display steps',ylegend='sea-level position')

We will assign 4 different colors based on relative sea-level change. To visualise the colors you can copy the html code below in the following website:

<!DOCTYPE html><html>
  <head>
    <style>
      #p1 {background-color:rgba(213,139,214,0.8);}
      #p2 {background-color:rgba(215,171,117,0.8);}
      #p3 {background-color:rgba(39,123,70,0.8);}
      #p4 {background-color:rgba(82,128,233,0.8);}
    </style>
  </head>
  <body>
    <p>RGB colors with opacity for Systems-Tracts:</p>
    <p id="p1">HST</p>
    <p id="p2">FSST</p>
    <p id="p3">LST</p>
    <p id="p4">TST</p>
  </body>
</html>

Define different colors for different systems tracts



In [ ]:

    
cHST = 'rgba(213,139,214,0.8)' 
cFSST = 'rgba(215,171,117,0.8)' 
cLST = 'rgba(39,123,70,0.8)' 
cTST = 'rgba(82,128,233,0.8)'

Specify the extent of every systems tracts according to the relative sea-level change



In [ ]:

    
HST1 = np.array([0,11],dtype=int)
FSST1 = np.array([11,26],dtype=int)
LST1 = np.array([26,30],dtype=int)
TST1 = np.array([30,36],dtype=int)
HST2 = np.array([36,42],dtype=int)
FSST2 = np.array([42,69],dtype=int)
LST2 = np.array([69,78],dtype=int)
TST2 = np.array([78,88],dtype=int)
HST3 = np.array([88,100],dtype=int)

Build the color list for each systems tract



In [ ]:

    
# Build the color list
STcolors = []
for k in range(HST1[0],HST1[1]):
    STcolors.append(cHST)
for k in range(FSST1[0],FSST1[1]):
    STcolors.append(cFSST)
for k in range(LST1[0],LST1[1]):
    STcolors.append(cLST)
for k in range(TST1[0],TST1[1]):
    STcolors.append(cTST)
for k in range(HST2[0],HST2[1]):
    STcolors.append(cHST)
for k in range(FSST2[0],FSST2[1]):
    STcolors.append(cFSST)
for k in range(LST2[0],LST2[1]):
    STcolors.append(cLST)
for k in range(TST2[0],TST2[1]):
    STcolors.append(cTST)
for k in range(HST3[0],HST3[1]):
    STcolors.append(cHST)

Plotting the resulting systems tracts



In [ ]:

    
# help(strata.viewSectionST)



In [ ]:

    
strata.viewSectionST(width = 1000, height = 500, cs=strat, colors=STcolors,
                   dnlay = 2, rangeX=[220000, 300000], rangeY=[-1050,-50], 
                   linesize=0.5, title='Systems tracts interpreted based on relative sea-level (RSL)')

2- Shoreline position / accomodation & sedimentation change

Before doing the interpretation using the two another methods, we need to calculate the shoreline positions (shoreline trajectory), the accommodation change ($\delta A$) and sedimentation change ($\delta S$) at shoreline position through time.
The endpoint of every deposited layer is also calculated to plot the Wheeler Diagram.
In order to calculate the shoreline trajectory, $\delta A$ and $\delta S$ at the shoreline, we need the information of the sedimentary layers when they deposited. In this example, 100 outputs are recorded (nbout). We use a loop to read the output every two timesteps (nstep).



In [ ]:

    
nbout = 100 
nstep = 2

The sed.timeT.pX.hdf5 file starts from time1. Therefore, the first file we load is T=1, then followed by T=2, 4, 6, 8, ..., 100



In [ ]:

    
strat_all = [strata.stratalSection(folder)]
strat_all[0].loadStratigraphy(1)
k = 1
for i in range(2,nbout+1,nstep):
    strat = strata.stratalSection(folder)
    strat_all.append(strat)
    strat_all[k].loadStratigraphy(i)
    k += 1

Building cross-sections parameters

We then slice the stratigraphic mesh to visualise the sedimentary architecture along a given cross-section.

To create the cross-section you will need to provide the position of the segment in the simulation space (xo,yo) and (xm,ym), a gaussian filter value for smoothing (gfilt a value of 0 can be used for non-smoothing) and the resolution of the cross-section (based on a number of points: nbpts).



In [ ]:

    
# Extract the value of sea-level for each stratigraphic layer
Time = np.array(time[::nstep])
sealevel = np.array(Sealevel[::nstep])
sealevel[0] = Sealevel[1]



In [ ]:

    
strata.viewData(x0 = Time/timeStep, y0 = sealevel, width = 600, height = 400, linesize = 3, 
                color = '#6666FF',title='Sea-level curve',xlegend='display steps',ylegend='sea-level position')



In [ ]:

    
npts = len(strat_all)
for i in range(0,npts):
    strat_all[i].buildSection(xo = x1, yo = y1, xm = x2, ym = y2, pts = nbpts, gfilter = gfilt)



In [ ]:

    
strat.buildParameters(npts,strat_all,sealevel)

3 - Shoreline trajectory method

Here we use a second method to interpret the depositional units based on the shifts of shoreline (from Helland-Hansen & Hampson (2009)):

Transgressive trajectory class TTC , when shoreline landward shift
Descending regressive trajectory class DRTC , when shoreline basinward shift accompanied by degradation
Ascending regressive trajectory class ARTC , when shoreline basinward shift accompanied by aggradation
Stationary trajectory class STC, when the shoreline stands still

Visualize shoreline trajectory and its gradient



In [ ]:

    
xval = Time/timeStep
yval = strat_all[0].dist[strat.shoreID_gs.astype(int)]
gval = np.gradient(yval)

# View shoreline position through time
#strata.viewData(x0 = xval, y0 = yval, width = 600, height = 400, linesize = 3, color = '#6666FF',
#               title='Shoreline trajectory',xlegend='display steps',ylegend='shoreline position in metres')

# Define the gradient evolution
#strata.viewData(x0 = xval, y0 = gval, width = 600, height = 400, linesize = 3, color = '#f4a142',
#               title='Shoreline trajectory gradient',xlegend='display steps',ylegend='gradient shoreline position')

We will assign 4 different colors based on relative shoreline change. To visualise the colors you can copy the html code below in the following website:

<!DOCTYPE html><html>
  <head>
    <style>
      #p1 {background-color:rgba(56,110,164,0.8);}
      #p2 {background-color:rgba(60,165,67,0.8);}
      #p3 {background-color:rgba(112,54,127,0.8);}
      #p4 {background-color:rgba(246,146,80,0.8);}
    </style>
  </head>
  <body>
    <p>RGB colors with opacity for Systems-Tracts:</p>
    <p id="p1">TC</p>
    <p id="p2">DRC</p>
    <p id="p3">ARC</p>
    <p id="p4">STC</p>
  </body>
</html>

Define different colors for different shoreline trajectory classes



In [ ]:

    
# Default color used: 
TC = 'rgba(56,110,164,0.8)'
DRC = 'rgba(60,165,67,0.8)'  
ARC = 'rgba(112,54,127,0.8)'
STC = 'rgba(252,149,7,0.8)'

# You can change them by specifying new values in the function
STcolors_ST = strata.build_shoreTrajectory(xval,yval,gval,sealevel,nbout,
                                           cTC=TC,cDRC=DRC,cARC=ARC,cSTC=STC)

Plotting the resulting shoreline trajectory classes



In [ ]:

    
strata.viewSectionST(width = 1000, height = 500, cs = strat, colors = STcolors_ST, 
                     dnlay = 2, rangeX=[220000, 300000], rangeY=[-1050,-50], 
                     linesize = 0.5, title = 'Classes interpreted based on shoreline trajectory (ST)')

4 - Accommodation succession method

This accommodation succession method is referred to Neal & Abreu (2009)), which is based on the relationship between the rate of changes in accommodation and sedimentation.

Retrogradation sequence set RSS
Aggradation to Progradation to Degradation (maybe) sequence set APDSS
Progradation to Aggradation sequence set PASS



In [ ]:

    
xval = Time/timeStep
ASval = strat.accom_shore-strat.sed_shore
gval = np.gradient(ASval)

# Accommodation (A) and sedimentation (S) change differences
#strata.viewData(x0 = xval, y0 = ASval, width = 600, height = 400, linesize = 3, 
#                color = '#6666FF',title='Change between accomodation & sedimentation',xlegend='display steps',
#                ylegend='A-S')

# Define the gradient evolution
#strata.viewData(x0 = xval, y0 = gval, width = 600, height = 400, linesize = 3, color = '#f4a142',
#               title='A&S gradient',xlegend='display steps',ylegend='gradient A&S')

We will assign 3 different colors based on relative accommodation space and sedimentation rate. To visualise the colors you can copy the html code below in the following website:

<!DOCTYPE html><html>
  <head>
    <style>
      #p1 {background-color:rgba(51,79,217,0.8);}
      #p2 {background-color:rgba(252,149,7,0.8);}
      #p3 {background-color:rgba(15,112,2,0.8);}
    </style>
  </head>
  <body>
    <p>RGB colors with opacity for Systems-Tracts:</p>
    <p id="p1">RSS</p>
    <p id="p2">PASS</p>
    <p id="p3">APDSS</p>
  </body>
</html>

Define different colors for different accommodation succession sequence sets.



In [ ]:

    
# Default color used: 
R = 'rgba(51,79,217,0.8)' 
PA = 'rgba(252,149,7,0.8)' 
APD= 'rgba(15,112,2,0.8)'

# You can change them by specifying new values in the function
STcolors_AS = strata.build_accomSuccession(xval,ASval,gval,nbout,cR=R,cPA=PA,cAPD=APD)

Plotting the resulting accommodation succession sequence sets



In [ ]:

    
strata.viewSectionST(width = 1000, height = 500, cs = strat, colors = STcolors_AS, dnlay = 2, 
                     rangeX=[220000, 300000], rangeY=[-1050,-50], linesize = 0.5, 
                     title = 'Sequence sets interpreted based on change of accommodation and sedimentation (AS)')

5 - Wheeler Diagram

Wheeler diagram (or chronostratigraphic chart) is a powerful tool to document unconformities between sequences, and to understand the stacking patterns of sedimentary successions and their relationship to sea level. It displays the horizontal distribution of contemporaneous sedimentary layer sequences, as well as hiatuses in sedimentation.

We define paleo-depositional environments (depositional facies) based on deposition depth. For example,

alluvial plain: 30m above sea level
shoreface: 0 - 50m below sea level
shallow marine: 50 - 200m below sea level
slope: 200 - 400m below sea level
deep marine: 400 - 600m below sea level
ocean basin: > 600m below sea level



In [ ]:

    
# Specify the depositional environment based on water depth
envIDs = [-30, 0, 50, 200, 400, 600]  
# Change the unit to be Myr
timeMA = [x/1e6 for x in Time] 
# For each depositional environments define a color 
color = ['limegreen','darkkhaki','sandybrown','khaki','c', 'teal', 'white']



In [ ]:

    
# help(strata.depthID)
# help(strata.viewWheeler)

Visualise the wheeler diagram



In [ ]:

    
# Get the position of depositional environments on each extracted layer
waterIDs = np.zeros((len(envIDs), npts))
for i in range(0,npts):
    waterIDs[:,i] = strata.depthID(cs = strat_all[i], sealevel = sealevel[i], envIDs = envIDs)
    
# Add the endpoint of deposition to the depositional environments extent
envi = np.append(waterIDs, [strat.depoend_gs], axis=0)

strata.viewWheeler(width = 10, height = 6, time = timeMA, rangeE = envi, shoreID = strat.shoreID, 
                dnlay = nstep, npts = npts, color = color, rangeX = [200, 300], rangeY = [-0.1,10.2], 
                linesize = 3, title = 'Wheeler Diagram', xlegend = 'Distance (km)', ylegend = 'Time (Myr)')

6 - Vertical stacking pattern

Vertical stacking patterns (synthetic wells) could also be extracted at different positions to show the thickness of parasequences and the trend of depositional units through time.

To extract the vertical staking, first we need to calculate the stacked strata thickness at each position, then color the strata according to the water depth when they are deposited.



In [ ]:

    
# help(strata.getStack)



In [ ]:

    
# position of wells /km
posit = [230, 240, 250, 260, 270]

colorFill, layTh = strata.getStack(cs = strat, posit = posit, envIDs = envIDs, color = color, dn = nstep)

Visualise the synthetic wells



In [ ]:

    
strata.viewStack(width = 4, height = 5, layTh = layTh, colorFill = colorFill)



In [ ]: