## Introduction

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 and sediment thickness , as well as the same location at some previous geological time when no sediments were present and the water depth was . The difference in water depth could have been caused by several reasons: the crustal deflection due to the sediment load , the global sea level change , and finally the thermo-tectonic subsidence .

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

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

From all the symbols in this equation the model expects , and (present-day sediment thickness) as inputs and it computes , and 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 , we reduce the pressure contribution from the mantle w.r.t. this level by . 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 ( is cancelled in the following equation):

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

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

## Sediment decompaction

Let’s now look more closely at how to compute sediment properties and . Bulk density is by definition

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

where 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 over time. Instead, we can get away with calculating only once from the present-day rock properties. But we still have to to find and for this we need a decompaction model.

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

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

Here 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 we can then compute at another burial depth by finding the zero of the following non-linear function:

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* function exhibits a smooth gradual and monotonic change over depth and hence a second-order Gauss quadrature is sufficient to approximate the integral:

The second required numerical technique is the Newton algorithm for solving :

## 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, , where is a surface porosity and 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,' 'maximum_burial']) 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'] @cached_property 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): self.events.update(events) def reconstruct_burial_history(self): """ Compute maximum burial depths for all the deposited layers. """ current_burial = 0.0 for event in self.events: 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 self.events[starting_event_id:]: thickness = event.layer.thickness_at_depth(current_burial) weight = event.layer.sediment_weight(constants) current_burial += thickness thickness_list.append(thickness) weight_list.append(weight) return thickness_list, weight_list def sea_level_change(self, event_id): """ Sea level change for a given event ID. """ return self.events[event_id].sea_level - self.initial_sea_level def bathymetry(self, event_id): """ Water depth value for a given event ID. """ return self.events[event_id].bathymetry

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(event_manager.events)): 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) subsidence.append(s) thickness_evolution.append(thickness) 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.