My current project involves developing a software tool that conducts the so-called backstripping analysis. In this post I would like to distill the information I’ve learned about this particular technique.

Backstripping simulation is a backward-in-time process where at each time step one removes the top-most (youngest) sediment layer together with its effects on the underlying sediments (via decompaction) and the Earth’s crust (isostasy). The purpose of this exercise is to compute the temporal evolution of the component of crustal subsidence that is not attributed to sediment loading. This subsidence is called thermo-tectonic, which essentially means that we have no clue as to why it happened but we suspect that it was caused by some combination of tectonic forces and thermal effects.

The idea behind backstripping is based on the following simplified viewpoint of the geological evolution. Let’s consider an abstract geological column at some location with water depth w and sediment thickness S, as well as the same location at some previous geological time when no sediments were present and the water depth was w_0. The difference in water depth could have been caused by several reasons: the crustal deflection due to the sediment load \xi, the global sea level change \Delta, and finally the thermo-tectonic subsidence T.


As you can see from my impeccably drawn cartoon, these variables can be linked together via

T = w + S - \xi - \Delta - w_0.

The w_0 term (which is omitted in other materials on backstripping) makes sure that T_0 = 0, which I find more intuitive: the model calculates incremental subsidence with respect to t = 0, so it makes sense to me to define T_0 as the zero baseline.

From all the symbols in this equation the model expects w(t), \Delta(t) and S_p (present-day sediment thickness) as inputs and it computes \xi(t), S(t) and T(t) as outputs.

Crustal deflection

For a 1D crustal deflection the simplest model is that of a hydrostatic equilibrium (or Airy isostasy). It says that there is a level deep down in the Earth’s mantle at which the hydrostatic pressure remains constant regardless of what happens at the surface. When we push the crust down by \xi, we reduce the pressure contribution from the mantle w.r.t. this level by \rho_mg\xi. According to the isostasy theory that same amount of pressure now needs to come from the load on top of the crust. We can see from the figure that the incremental crustal load has three components: sediment load, water load due to deflection and water load due to sea level change (g is cancelled in the following equation):

\rho_m\xi = \rho_bS + \rho_w(\xi-S) + \rho_w\Delta,

\xi = \frac{\rho_b-\rho_w}{\rho_m-\rho_w}S + \frac{\rho_w}{\rho_m-\rho_w}\Delta

By the way, \rho_m, \rho_b and \rho_w are densities of mantle, sediment column (b is for bulk density) and water respectively. When we substitute this expression for \xi back into the backstripping equation, we obtain the final result:

T = w - w_0 + \frac{\rho_m - \rho_b}{\rho_m - \rho_w}S - \frac{\rho_m}{\rho_m - \rho_w}\Delta.

Note that these equations also apply when sediments are exposed to air instead of water (e.g. mountains), if we set \rho_w = 0.

Sediment decompaction

Let’s now look more closely at how to compute sediment properties S(t) and \rho_b(t). Bulk density is by definition

\rho_b(t) = \phi(t)\rho_w + [1-\phi(t)]\rho_s,

where \phi is porosity and \rho_s is sediment rock density (e.g. pure quartz). Notice, that the corresponding term in the expression for \xi is

(\rho_b(t)-\rho_w)S(t) = (\rho_s - \rho_w)[1-\phi(t)]S(t) = (\rho_s - \rho_w)h,

where h = [1-\phi(t)]S(t) is the “rock thickness” of the sediment column which does not change over time if we neglect chemical reactions. That is, a sediment layer may compact over time due to grains being squeezed together from the weight above, but the actual mass of rock remains the same.

This is good news because it means that there is no need to calculate \rho_b(t) over time. Instead, we can get away with calculating h only once from the present-day rock properties. But we still have to to find S(t) and for this we need a decompaction model.

To derive it we will rely again on the fact that h is constant over time. We will also assume that the temporal change of porosity can be captured with a porosity-depth relationship \phi(z). In other words, porosity changes over time only due to the sediments progressive burial and hence \phi(t) = \phi(z(t)).

Let’s consider an arbitrary sediment layer of thickness S_p and calculate its “rock thickness” h from the present-day stratigraphic information:

h = \int_{z_{\mathrm{max}}}^{z_{\mathrm{max}}+S_p}[1-\phi(z)]\mathrm{d}z = S_p - \int_{z_{\mathrm{max}}}^{z_{\mathrm{max}}+S_p}\phi(z)\mathrm{d}z.

Here z_{\mathrm{max}} is the maximum burial depth of our layer, which does not necessarily equal its present-day burial depth if a major erosional event caused its uplift sometime in the geologic past. From knowing h we can then compute S(z) at another burial depth z by finding the zero of the following non-linear function:

f(S) = S - \int_{z}^{z+S}\phi(z)\mathrm{d}z - h,

which enforces sediment mass concentration during compaction.

In practical terms, we need to employ two numerical methods to model sediment compaction. First is the numerical integration of a porosity-depth relationship. A typical theoretical \phi(z) function exhibits a smooth gradual and monotonic change over depth and hence a second-order Gauss quadrature is sufficient to approximate the integral:

\Delta z = \frac{z_2-z_1}{2}, \quad \bar{z} = \frac{z_1+z_2}{2},

\int_{z_1}^{z_2}\phi(z)\mathrm{d}z\approx\Delta z\left[\phi\left(\bar{z}-\frac{\sqrt{3}}{3}\Delta z\right)+\phi\left(\bar{z}+\frac{\sqrt{3}}{3}\Delta z\right)\right].

The second required numerical technique is the Newton algorithm for solving f(S)=0:

S^{(0)}=h,\quad S^{(j+1)}=S^{(j)}-\frac{f(S^{(j)})}{f^{'}(S^{(j)})}=S^{(j)}-\frac{f(S^{(j)})}{1-\phi(z_{i}+S^{(j)})}.

Sample implementation

Here I want to present a sample Python implementation of the described model in case anyone finds it useful. For our implementation we will use Athy’s porosity-depth relationship, \phi(z) = \phi_0 \exp(-cz), where \phi_0 is a surface porosity and c is a compaction rate. We will only consider depositional events, that is, once a layer had been deposited, it remained in place from then on without being eroded.

We first implement the decompaction related functionality in the Layer class:

import math

from traits.api import (HasStrictTraits, Function, Float,
                        Property, cached_property)

class Layer(HasStrictTraits):
    maximum_burial = Float
    present_thickness = Float
    porosity_function = Function
    compaction_rate = Float
    surface_porosity = Float
    sediment_density = Float
    sediment_thickness = Property(depends_on=['present_thickness,'

    def set_rock_properties(self, properties_dict):
        self.surface_porosity = properties_dict['surface_porosity']
        self.compaction_rate = properties_dict['compaction_rate']
        self.sediment_density = properties_dict['sediment_density']

    def _get_sediment_thickness(self):
        """ Compute sediment thickness. """
        z0 = self.maximum_burial
        z1 = z0 + self.present_thickness
        water_thickness = self.integrate_porosity_function(z0, z1)
        return self.present_thickness - water_thickness

    def integrate_porosity_function(self, z0, z1):
        """ Numerically integrate porosity function over the given interval. """
        w = 0.5773502691896257  # sqrt(3)/3
        halflength = 0.5 * (z1 - z0)
        midpoint = 0.5 * (z0 + z1)

        porosity_0 = self.porosity_function(self, midpoint + halflength * w)
        porosity_1 = self.porosity_function(self, midpoint - halflength * w)
        return halflength * (porosity_0 + porosity_1)

    def thickness_at_depth(self, depth, eps=1e-6):
        """ Computes layer's thickness if buried at a given depth. """
        thickness = self.present_thickness  # initial guess
        # Newton iteration
        carry_on = True
        while carry_on:
            water_thickness = self.integrate_porosity_function(depth, depth + thickness)
            function_value = thickness - self.sediment_thickness - water_thickness
            derivative_value = 1.0 - self.porosity_function(self, depth + thickness)
            thickness -= function_value / derivative_value
            carry_on = abs(function_value) > eps
        return thickness

    def sediment_weight(self, constants):
        """ Layer weight above that of water. """
        return (self.sediment_density - constants.water_density) \
            * constants.gravity * self.sediment_thickness

def athy_porosity(layer, z):
    """ Athy's porosity-depth relationship. """
    return layer.surface_porosity * math.exp(-layer.compaction_rate * z)

The Layer objects know how to compute their rock thickness and how thick they would be at different burial depths. To describe a collection of layers deposited at some time with a certain geologic context, we also introduce the Deposition class, as well as the EventManager class which holds a collection of deposition events and can perform calculations which require knowledge of the burial history:

from sortedcontainers import SortedListWithKey
from traits.api import HasStrictTraits, Float, Instance

from layer import Layer

class Deposition(HasStrictTraits):
    age = Float
    bathymetry = Float
    sea_level = Float
    layer = Instance(Layer)

class EventManager(HasStrictTraits):
    events = Instance(SortedListWithKey, kw={'key': lambda e: e.age})
    initial_age = Float
    initial_sea_level = Float
    initial_bathymetry = Float

    def add_events(self, events):

    def reconstruct_burial_history(self):
        """ Compute maximum burial depths for all the deposited layers. """
        current_burial = 0.0
        for event in
            event.layer.maximum_burial = current_burial
            current_burial += event.layer.present_thickness

    def decompact_layers(self, starting_event_id, constants):
        """ Decompaction of a sediment column. """
        current_burial = 0.0
        thickness_list = []
        weight_list = []
        for event in[starting_event_id:]:
            thickness = event.layer.thickness_at_depth(current_burial)
            weight = event.layer.sediment_weight(constants)
            current_burial += thickness

        return thickness_list, weight_list

    def sea_level_change(self, event_id):
        """ Sea level change for a given event ID. """
        return[event_id].sea_level - self.initial_sea_level

    def bathymetry(self, event_id):
        """ Water depth value for a given event ID. """

Assuming that we have a setup EventManager object, the actual backstripping algorithm is trivial to implement:

def compute_deflection(sediment_weight, sea_level_change, constants):
    """ helper function for Airy isostasy. """
    total_weight = sediment_weight + constants.gravity * constants.water_density * sea_level_change
    return total_weight / (constants.gravity * (constants.mantle_density - constants.water_density))

def compute_subsidence(event_manager, constants=PhysicalConstants()):
    """ Actual backstripping is performed here. """
    subsidence = []
    thickness_evolution = []
    for event_id in range(len(
        thickness, weight = event_manager.decompact_layers(event_id, constants)
        total_thickness = sum(thickness)
        total_weight = sum(weight)
        sea_level_change = event_manager.sea_level_change(event_id)
        bathymetry = event_manager.bathymetry(event_id)
        deflection = compute_deflection(total_weight, sea_level_change, constants)
        s = (bathymetry + total_thickness - deflection - sea_level_change
             - event_manager.initial_bathymetry)
    return subsidence[::-1], thickness_evolution[::-1]

The compute_subsidence function returns the subsidence values as well as the whole sediment column thickness distribution for each geologic age, sorted from earliest to present-day.

Example application

The example application is taken from the Basin Analysis book by Allen and Allen (2013). Here is the main file listing that contains all the data:

from backstripping import prepare_events, compute_subsidence, plot_results

# setup layers
rock_properties = [
    {'surface_porosity': 0.63, 'compaction_rate': 0.51e-3, 'sediment_density': 2720.0},
    {'surface_porosity': 0.49, 'compaction_rate': 0.27e-3, 'sediment_density': 2650.0},
    {'surface_porosity': 0.70, 'compaction_rate': 0.71e-3, 'sediment_density': 2710.0},
    {'surface_porosity': 0.40, 'compaction_rate': 0.60e-3, 'sediment_density': 2720.0},
    {'surface_porosity': 0.20, 'compaction_rate': 0.60e-3, 'sediment_density': 2870.0},
    {'surface_porosity': 0.05, 'compaction_rate': 0.20e-3, 'sediment_density': 2960.0},

ages = [260, 245, 210, 160, 145, 125, 100, 80, 55, 45, 0]
sea_levels = [10, 0, 0, -20, -40, 70, 80, 100, 50, 40, 0]
bathymetries = [-20, 0, 20, 10, 20, 20, 200, 300, 350, 325, 300]
rock_types = [4, 5, 1, 4, 3, 1, 2, 0, 1, 0]
thicknesses = [400, 750, 250, 400, 200, 900, 1300, 750, 250, 200]

event_manager = prepare_events(ages, bathymetries, sea_levels, thicknesses, rock_types, rock_properties)
subsidence, thickness_evolution = compute_subsidence(event_manager)
plot_results(ages, subsidence, thickness_evolution, sea_levels, bathymetries)

And here is the result:


As you can see, less than half of the present-date crust depression for this location is due to tectonic forces and most of it is instead caused by the weight of deposited sediments.

All the code necessary to run this example and generate this plot can be found here.