PALEOSTRIPv1.0 – a user-friendly 3D backtracking software to reconstruct paleo-bathymetries

Paleo-bathymetric reconstructions provide boundary conditions to numerical models of ice sheet evolution and ocean circulation that are critical to understanding their evolution through time. The geological community lacks a complex open-source tool that allows for community implementations and strengthens research synergies. To fill this gap, we present PALEOSTRIPv1.0, a MATLAB open-source software designed to perform 1D, 2D, and 3D backtracking of paleo-bathymetries. PALEOSTRIP comes with a graphical user interface (GUI) to facilitate computation of sensitivity tests and to allow the users to switch all the different processes on and off and thus separate the various aspects of backtracking. As such, all physical parameters can be modified from the GUI. It includes 3D flexural isostasy, 1D thermal subsidence, and possibilities to correct for prescribed sea level and dynamical topography changes. In the following, we detail the physics embedded within PALEOSTRIP, and we show its application using a drilling site (1D), a transect (2D), and a map (3D), taking the Ross Sea (Antarctica) as a case study. PALEOSTRIP has been designed to be modular and to allow users to insert their own implementations.

1 Introduction

Ongoing climate changes are urging the scientific community to project future climate evolution in response to carbon emission trajectories (e.g. the Shared Socio-economical Pathways by Riahi et al., 2017). The Coupled Model Intercomparison Project (CMIP), now ending phase 6 (Eyring et al., 2016), has been producing a large amount of climate projections that extend to 2100. However, some of the climatic variables, e.g. the deep ocean, the carbon cycle, and the ice sheets and glaciers, react more slowly to climate changes (e.g. Colleoni et al., 2018b; Noble et al., 2020), despite already showing evidence of changes over the past few decades (Caesar et al., 2021). Their main response has yet to be observed and is likely to happen beyond the 21st century, which is encouraging climatologists to project changes on longer timescales that cover centuries to millennia into the future (e.g. IPCC, 2013, 2019). Projecting to such timescales implies designing corresponding realistic carbon emission trajectories, and so far millennial-scale emissions trajectories are just extension of existing emission scenarios beyond 2100 (e.g. Golledge et al., 2015; DeConto and Pollard, 2016). At this point, the exercise becomes difficult, and this is when reconstructing past climates becomes important. It allows for the testing of the response of the Earth's climate under different but realistic atmospheric greenhouse gas (GHG) concentrations (e.g. Haywood et al., 2019). GHG concentrations higher than present-day or future levels, i.e. larger than 400 ppm for atmospheric CO₂, can only be found for times prior to 3 million years ago (e.g. Zhang et al., 2013; Bracegirdle et al., 2019). Going back to those times (or even before), it is likely that the tectonic setting responsible for other boundary conditions, such as oceanic gateways, elevation of mountain ranges, continental margin expansion, and the location and extent of continental masses themselves, differed from that of the present-day (e.g. Herold et al., 2008; Frigola et al., 2018; Müller et al., 2018a; Straume et al., 2020; Hochmuth et al., 2020). Nevertheless, simulating and reconstructing past climatic conditions can bring useful hints as to how the future climate might evolve and also help in narrowing the range of likely long-term carbon emissions trajectories.

Reconstructing past topographies and bathymetries is fundamental for paleoclimate and paleo-ice-sheet simulations (e.g. Otto-Bliesner et al., 2017; Colleoni et al., 2018a; Straume et al., 2020; Paxman et al., 2020). The numerous oceanic deep-drilling campaigns that occurred over the past decades have the potential to constrain such reconstructions, but this is not the case when part of the information has been eroded and/or reworked during geological time by other tectonic or climatic processes. In addition, during sediment deposition, the bathymetry itself changes due to different processes, such as the loading of accumulated sediments, thermal subsidence acting on extended continental crust, or subsidence or uplift induced by mantle dynamics (dynamic topography) or sea level changes (e.g. Kirschner et al., 2010; Celerier, 1988; Sclater and Christie, 1980). Other external factors also influence the bathymetry, e.g. ice loading in polar areas. When accounting for all those factors and processes by decompacting and removing overlying sediments, it is possible to backtrack the bathymetry, and thus the paleo-water depth, of a chosen specific time interval. Conversely, if the target of the study is to reconstruct the burial history of a sedimentary basin, the technique is the same but needs constraints on the paleo-water depths and is called “backstripping” (Steckler and Watts, 1978).

There are few existing open-source backstripping or backtracking codes. Flex-Decomp by Badley Geoscience Ltd. (Kusznir et al., 1995) allows for 2D-flexural backstripping or backtracking, but its code is not open-source. A 3D version of Flex-Decomp exists, but it is not open-source (Roberts et al., 2003) and is used exclusively by Badley Geoscience Ltd. and is thus not available externally. It comes along with the sister programme Stretch by Badley Geoscience Ltd., a software used to compute forward modelling of basin evolution that provides a spatially variable stretching factor cross section to be used by Flex-Decomp. BasinVis (Lee et al., 2016, 2020) is a MATLAB open-source code with a graphical user interface. It calculates compaction trends from input drill site (1D) lithological units. The estimated compaction trends can be applied in thickness restoration of stratigraphic units and subsidence analysis data that can be spatially interpolated between input drill sites to reconstruct a temporal basin evolution. Dressel et al. (2015, 2017) perform 3D backtracking of the southwestern African continental margin while accounting for 1D thermal subsidence and Pratt isostasy, taking into account lateral density variations. PyBacktrack is a 1D backtracking and backstripping open-source code (Müller et al., 2018b) aimed at reconstructing paleo-bathymetries. It allows the processing of drilling sites both on oceanic and continental crust; can be connected to the suite of geodynamical open-source software GPlates (https://www.gplates.org/, last access: August 2021); and benefits from geodynamical corrections related to kinematic, tectonic, and geodynamic models of tectonic plate movements through time. DeCompactionTool (Hölzel et al., 2008) proposes a similar approach to pybacktrack (also in 1D) but allows for the performing of a Monte Carlo style analysis, i.e. performing a large number of 1D runs based on a possible range of main physical parameters defined by admissible minimum and maximum values to provide a quantitative estimate of the backstripping error.

Both 3D flexural backstripping and backtracking are needed to reconstruct basin-wide or continental-wide areas that will be prescribed as boundary conditions within climate and ice sheet models. A few studies mentioned the use of 3D backstripping (e.g. Scheck et al., 2003; Hansen et al., 2007; Smallwood, 2009; Amadori et al., 2017), and a very limited number of studies provide this with 3D flexural backstripping (Roberts et al., 2003; Steinberg et al., 2018; Roberts et al., 2019).

Here we present PALEOSTRIP, a MATLAB open-source software designed to perform 1D, 2D, and 3D backtracking of paleo-bathymetries. PALEOSTRIP comes with a graphical user interface (GUI) to facilitate computation of sensitivity tests and allows the users to switch all of the different processes on and off and thus separate the various aspects of backtracking. As such, all physical parameters can be modified from the GUI. It includes 3D flexural isostasy, 1D thermal subsidence, and possibilities to correct for prescribed sea level and dynamical topography changes. In the following, we detail the physics embedded within PALEOSTRIP and show a few applications for a drilling site (1D), a transect (2D), and a map (3D), taking the Ross Sea (Antarctica) as a case study.

Figure 1 Example of the PALEOSTRIP application GUI, showing this following elements: (a) the main PALEOSTRIP GUI, (b) tab interface to set thermal subsidence (bottom of main GUI), (c) tab interface to set sea level corrections, and (d) tab interface to set dynamic topography correction.

Figure 2 The backtracking procedure used in PALEOSTRIP is as follows. Sediment layers L are separated by well-defined seismic or lithological unconformities U. All the sediment layers have different lithological properties defined by their porosity ϕ and density ρ. These two properties change during the decompaction process, departing from the present-day depth of sediment layers, as buried sediments are unloaded from the upper sediment layers (isostatic correction). The water depth at each time step is calculated using a time-varying thermal subsidence model (Eq. 3).

2 Model framework and requirements

PALEOSTRIP is a MATLAB open-source code developed under the GNU General Public License v3.0. It is composed of a set of routines accessed using a graphical interface from which users can load data, change physical parameters, and plot and save results (Fig. 1). The code is distributed on GitHub (https://github.com/flocolleoni/PALEOSTRIPv1.0, last access: August 2021) alongside a user manual providing all necessary explanations and examples of how to input data and use PALEOSTRIP functionalities. Note that the final version related to this study has been corrected for minor bugs and that the main “plot and save” graphical interface has been modified to implement 12 colour scales, including 9 colour-blind-friendly scales from Crameri et al. (2020). We thus again invite the readers to download the code from GitHub and from Zenodo.

PALEOSTRIP has been designed and coded with MATLAB (2019) and runs on any operating system. Most of the code should be compatible with previous versions of MATLAB, provided that it includes the App Designer toolbox (https://it.mathworks.com/products/matlab/app-designer.html, last access: August 2021). The code is incompatible with GUIDE and from MATLAB R2020 onwards GUIDE is no longer available. PALEOSTRIP cannot be run without its graphical interface, and thus a version of MATLAB with the App Designer toolbox is necessary.

2.1 PALEOSTRIP graphical user interface

PALEOSTRIP graphical user interface (GUI) is composed of two different tabs. The first one is dedicated to input data, physical parameters, and choices of physical methods for each of the processes accounted for in backtracking (Fig. 1). The second one is dedicated to plotting backtracked data and saving results (see Sect. 6). The GUI can be launched either through the MATLAB main interface and double-clicking on the corresponding App Designer GUI file or by exporting it as a stand-alone application that can be executed without opening MATLAB. We do not provide a built-in PALEOSTRIP application, since the compilation of the applications depends on the operating system on which it is created, but we provide the PALEOSTRIP code instead. As such, the user is free to export it from the App Designer tool.

2.2 Format of input and output data

2.2.1 Coordinate system

Cartesian coordinates (in m) are required to perform flexural isostasy calculations, as the flexural response depends on the distance from the load. Consequently, input data must be provided on a Cartesian coordinate grid. PALEOSTRIP does not run on geographical coordinates and does not provide any tool to convert input data from geographical to Cartesian coordinates. However, many examples of open-source software or code exist to perform this step, such as the MATLAB Mapping Toolbox (https://it.mathworks.com/products/mapping.html, last access: August 2021), Generic Mapping Tools 6 (https://www.generic-mapping-tools.org/, last access: August 2021), OBLIMAP2 (Reerink et al., 2016), or any other geographic information system (GIS) software.

2.2.2 Input files

Horizon depths have to be provided individually in single ASCII files, defining the given quantity along with its coordinates. For example, horizon depth Z for a 1D drill site, extension along the transect and depth (X,Z) for 2D transects, and the horizontal position (X,Y) and depth (Z) for 3D horizon maps. Output data are saved with the same format. PALEOSTRIP does not read and write grids in the NetCDF format. Lithological parameters must be provided in a separate ASCII file and are spatially uniform. In the present study, the input data files associated with each case study are zipped in paleostrip_examples.zip, available on GitHub at (https://github.com/flocolleoni/PALEOSTRIPv1.0, last access: August 2021). Since the lithological parameters vary substantially given the composition of sediments and their depositional context, we refer the reader to Kominz et al. (2011) for values and detailed explanations on how to retrieve the decompaction coefficients for the different lithologies of marine sediments.

Figure 3 Workflow of the code implemented in PALEOSTRIP illustrating the various steps of computation according to selected options in the GUI. N sediment layers are decompacted during the backtracking procedure. Smooth blue rectangles correspond to the start and end of the backtracking run. Light red parallelograms indicate input and output data. Light green rectangles are intermediate steps in the backtracking computation. Green diamonds correspond to options selected by the user in the GUI. If isostasy, thermal subsidence, sea level changes, or dynamic topography are switched off, the corresponding correction equals zero.

Figure 4 The 2D and 3D domain expansions for computation of flexural isostasy. (a) In the case of 2D vertical sections, input data are extended with the edge values for 90 % of the initial total array length (NX). This extended array is placed at the centre of a null array that is 3 times larger than the initial input data array. (b) In the case of 3D maps, input data are extrapolated with a nearest neighbour algorithm on a mesh grid (black) to obtain a structured squared grid if the initial data are given on an irregular polygon. This mesh grid is extended by 30 % on all edges (red) and is placed at the centre of a grid twice as large as the mesh grid obtained from initial input data (blue).

3 PALEOSTRIP: backtracking and limitations

During its evolution, a submarine continental margin can experience various processes that modify its morphology. For example, a rifted basin can form in response to plate tectonics displacement and/or erosion that shapes the basic structure of the margin. Surface erosion from the hinterland can supply the margin with sediments that fill morphological depressions if the accommodation space is large enough. The accommodation space depends on initial conditions, on the tectonics, and on the thermal subsidence of the continental margin, i.e. the lithosphere and asthenosphere cooling during and after the rifting that leads to a deepening of the margin over time. Eustatic and regional sea level changes modulate the available accommodation space and the distance from the sediment sources. To reconstruct the past subsidence history of a continental margin at a given time, all those processes have to be accounted for and corrected in a procedure called “backstripping” (Steckler and Watts, 1978). Backstripping consists of decompacting and removing sediment layers iteratively back in time to reconstruct the past tectonic subsidence history of a basin or a margin (Fig. 2):

$$\begin{array}{}\text{(1)}& {\displaystyle }{\displaystyle }Z\left(i\right)={\displaystyle \frac{{\mathit{\rho }}_{\mathrm{m}}-{\mathit{\rho }}_{\mathrm{s}}}{{\mathit{\rho }}_{\mathrm{m}}-{\mathit{\rho }}_{\mathrm{w}}}}S\left(i\right)+\text{WD}\left(i\right)-{\displaystyle \frac{{\mathit{\rho }}_{\mathrm{m}}}{{\mathit{\rho }}_{\mathrm{m}}-{\mathit{\rho }}_{\mathrm{w}}}}\mathrm{\Delta }\text{SL}\left(i\right),\end{array}$$

where the first term of the right-hand side corresponds to the isostatic compensation (ρ_s sediment density, ρ_m mantle density, ρ_w seawater density) of the ith sediment layer thickness S(i) accumulated on the basement, WD(i) is the paleo-water depth of the ith layer, and the last term of the equation corresponds to the water load correction due to sea level variations (ΔSL) at time of the ith layer. This equation solves the time evolution of total subsidence Z(i) of the basement, and WD needs to be provided for each i layer to solve the equation. Conversely, in cases where the focus is on reconstructing time-varying paleo-water depths WD(i), the procedure is called “backtracking”, and a total subsidence Z(i) time evolution needs to be provided as input to the equation:

$$\begin{array}{}\text{(2)}& {\displaystyle }{\displaystyle }\text{WD}\left(i\right)=Z\left(i\right)-{\displaystyle \frac{{\mathit{\rho }}_{\mathrm{m}}-{\mathit{\rho }}_{\mathrm{s}}}{{\mathit{\rho }}_{\mathrm{m}}-{\mathit{\rho }}_{\mathrm{w}}}}S\left(i\right)+{\displaystyle \frac{{\mathit{\rho }}_{\mathrm{m}}}{{\mathit{\rho }}_{\mathrm{m}}+{\mathit{\rho }}_{\mathrm{w}}}}\mathrm{\Delta }\text{SL}\left(i\right).\end{array}$$

PALEOSTRIP is a backtracking software and is designed to reconstruct paleo-water depths given a provided subsidence history.

The equations above are the original backstripping and backtracking equations developed by Steckler and Watts (1978) and are explained therein in major detail. Over the past few decades, it has also been found that mantle convection generates changes in the regional topography, i.e. so-called “dynamic topography” (Müller et al., 2018a). Paleo-water depth needs to be corrected for dynamic topography changes ΔDynT(i) causing uplift or subsidence of the topography and bathymetry:

$$\begin{array}{ll}{\displaystyle }\text{WD}\left(i\right)=& {\displaystyle }Z\left(i\right)-{\displaystyle \frac{{\mathit{\rho }}_{\mathrm{m}}-{\mathit{\rho }}_{\mathrm{s}}}{{\mathit{\rho }}_{\mathrm{m}}-{\mathit{\rho }}_{\mathrm{w}}}}S\left(i\right)+{\displaystyle \frac{{\mathit{\rho }}_{\mathrm{m}}}{{\mathit{\rho }}_{\mathrm{m}}+{\mathit{\rho }}_{\mathrm{w}}}}\mathrm{\Delta }\text{SL}\left(i\right)\\ \text{(3)}& {\displaystyle }& {\displaystyle }+\mathrm{\Delta }\text{DynT}\left(i\right).\end{array}$$

In this equation, the second and third terms on the right-hand side have been written accounting for Airy local isostatic compensation. However, PALEOSTRIP also makes use of 2D and 3D flexural isostasy (see Sect. 3.2). In this case, the second and third terms on the right-hand side of the equation represent the depth correction due to the flexural response to unloading of sediment during backtracking and to water loading or unloading due to sea level changes.

In the following, we explain how each term of Eq. (3) is treated within PALEOSTRIP. Most of the equations reported below are taken from Allen and Allen (2013) if not otherwise specified. The various aspects of the backtracking procedure are explained as part of the PALEOSTRIP workflow (Fig. 3).

Note that the physics implemented do not allow for the treatment of oceanic crust in this version. This can be done by adding a few more options in the GUI, mainly for thermal subsidence (e.g. Müller et al., 2018b). This can be easily implemented.

3.1 Decompaction

The total sediment thickness S accumulated at a given time can be obtained by reconstructing its compaction history. Eroded sediments are transported to the margin deposit and accumulate wherever the accommodation space allows for it. Accumulated sediments compact over time under loading, e.g. by overlying sediments. Therefore, to reconstruct the paleo-architecture of a continental margin at a given time, the sediment layers of a margin are stripped off sequentially, and the remaining underlying sediments need to be gradually decompacted. Decompaction is thus central to backstripping and backtracking.

Compaction or decompaction of sediments both imply a change in total volume V_total of deposited sediments; this is mostly caused by changes in the porosity of sediments and (to a lesser extent) changes in sediment compression. In submarine environments, sediments are saturated by water, and their compaction implies a decrease in pore fluid pressure via expelling of the water out of the sediment layers. Conversely, decompaction involves an increase in the pore fluid pressure by injecting water within the compacted sediments. The decompaction process consists of calculating the changes in porosity of the various sediment layers to determine the amount of water that was contained in the sediments at the time of their deposition based on their lithology. The depth of the decompacted sediment layers is recalculated on the basis of those porosity changes. Laboratory experiments have determined that for large depths, the evolution of porosity can be described by an exponential relationship rather than by a linear empirical equation:

$$\begin{array}{}\text{(4)}& {\displaystyle }{\displaystyle }\mathit{\varphi }={\mathit{\varphi }}_{\mathrm{0}}{e}^{-cz},\end{array}$$

where ϕ₀ corresponds to deposition porosity at the seafloor (or at surface if emerged), c is the exponential slope decay coefficient (also referred to here as the compaction coefficient) depending on the lithology expressed in 1 km⁻¹, and z is depth in kilometres. This empirical equation implies that ϕ₀ is decreased by the factor of $\mathrm{1}/e$ at the depth of $\mathrm{1}/c$ km.

In general, the mass of the total sediment column at a given time does not change: in submarine environments the water that is expelled from the sediments is implicitly added to the water column above the underlying sediments during the compaction process. Thus, only the volume of the sediment layers V_total changes due to water volume changes V_w, which allows for the integration of the porosity Eq. (4). Considering a unit's cross-sectional area, the thickness of water W_T to be added to the compacted sediment volume in order to decompact it is given by the following equation:

$$\begin{array}{}\text{(5)}& {\displaystyle }{\displaystyle }{W}_{\mathrm{T}}=\underset{z_i}{\overset{z_{i+\mathrm{1}}}{\int }}{\mathit{\varphi }}_{\mathrm{0}}{e}^{-cz}\mathrm{d}z,\end{array}$$

which upon integration gives

$$\begin{array}{}\text{(6)}& {\displaystyle }{\displaystyle }{W}_{\mathrm{T}}\left(i\right)={\displaystyle \frac{{\mathit{\varphi }}_{\mathrm{0}}}{c}}\mathit{\{}{e}^{-c{z}_i^{\prime }}-e^{-c{z}_{i+\mathrm{1}}^{\prime }}\mathit{\}},\end{array}$$

where $z_i^{\prime }$ and $z_{i+\mathrm{1}}^{\prime }$ are the newly decompacted depths of z_i and z_i+1. Because the volume of sediment grains V_sed remains unchanged, the integrated compacted sediment thickness S_dt between two given depths, considering a unit's cross-sectional area, can be written as follows:

$$\begin{array}{}\text{(7)}& {\displaystyle }{\displaystyle }{S}_{\text{dt}}=z_{i+\mathrm{1}}-z_i-{\displaystyle \frac{{\mathit{\varphi }}_{\mathrm{0}}}{c}}\mathit{\{}{e}^{-c{z}_i}-e^{-c{z}_{i+\mathrm{1}}}\mathit{\}},\end{array}$$

Based on the fact that $V_{\text{total}}=V_{\text{sed}}+V_{\mathrm{w}}$, the decompacted depths z^′(i) for a unit's cross-sectional area can be inferred, and the final decompaction equation is given by the following equation:

$$\begin{array}{ll}{\displaystyle }{z}_{i+\mathrm{1}}^{\prime }-z_i^{\prime }=& {\displaystyle }{z}_{i+\mathrm{1}}-z_i-{\displaystyle \frac{{\mathit{\varphi }}_{\mathrm{0}}}{c}}\mathit{\{}{e}^{-c{z}_i}-e^{-c{z}_{i+\mathrm{1}}}\mathit{\}}\\ \text{(8)}& {\displaystyle }& {\displaystyle }+{\displaystyle \frac{{\mathit{\varphi }}_{\mathrm{0}}}{c}}\mathit{\{}{e}^{-c{z}_i^{\prime }}-e^{-c{z}_{i+\mathrm{1}}^{\prime }}\mathit{\}}.\end{array}$$

PALEOSTRIP solves this equation iteratively (starting with $z_i^{\prime }=\mathrm{0}$ and then the further iterations) until convergence is reached. Lithological parameters, i.e. the decompaction coefficient c and the depositional surface porosity ϕ₀, are prescribed in the user-provided parameter files (see Sect. 6).

3.2 Isostatic correction

Sediments load the basement of the continental margin and cause a local and regional subsidence during accumulation. Water also loads the basement due to eustatic sea level variations, leading to changes in the water depth WD that are also caused by ocean bottom subsidence and/or uplift through time. During the decompaction, the newly computed decompacted depths of the remaining sediment layers also need to be adjusted to account for the changes in their porosity and thus their density due to the presence of water filling the pores. For each sediment layer, the porosity ρ_ds(i) of the decompacted thickness S_d(i) of the ith sediment layer is

$$\begin{array}{}\text{(9)}& {\displaystyle }{\displaystyle }{\mathit{\varphi }}_{\text{ds}}\left(i\right)={\displaystyle \frac{{\mathit{\varphi }}_{\mathrm{0}}}{c}}{\displaystyle \frac{e^{-c{z}^{\prime }(i+\mathrm{1})}-e^{-c{z}^{\prime }\left(i\right)}}{S_{\mathrm{d}}\left(i\right)}},\end{array}$$

and the decompacted sediment bulk density ρ_ds(i) of the ith sediment layer is given by

$$\begin{array}{}\text{(10)}& {\displaystyle }{\displaystyle }{\mathit{\rho }}_{\text{ds}}\left(i\right)={\mathit{\varphi }}_{\text{ds}}\left(i\right){\mathit{\rho }}_{\mathrm{w}}+(\mathrm{1}-{\mathit{\varphi }}_{\text{ds}}(i\left)\right){\mathit{\rho }}_S\left(i\right).\end{array}$$

At each time step, the removal of the top layer causes the decompaction of the underlying sediment layers.

Removed sediment layers are substituted with seawater. Unloading causes an isostatic compensation that can be calculated by using different isostatic methods. In PALEOSTRIP, two methods are implemented.

3.2.1 Airy local compensation

The Airy local compensation is the most used isostatic method in backstripping and backtracking. It involves a local depth compensation Z_airy(i) for which the rocks and sediments and the underlying asthenosphere are considered to be in hydrostatic equilibrium:

$$\begin{array}{}\text{(11)}& {\displaystyle }{\displaystyle }{Z}_{\text{airy}}\left(i\right)=S_{\mathrm{d}}\left(i\right)\left({\displaystyle \frac{{\mathit{\rho }}_{\mathrm{m}}-{\mathit{\rho }}_{\text{ds}}\left(i\right)}{{\mathit{\rho }}_{\mathrm{m}}-{\mathit{\rho }}_{\mathrm{w}}}}\right).\end{array}$$

This method implies that the weight of the sediments in a given grid point does not impact the adjacent points. Therefore, in the case of a 2D transect or a 3D map, grid points are independent of each other. This is of course not realistic, and the Airy local compensation preferably should be applied to the decompaction of wells only (Roberts et al., 1998).

3.2.2 Flexural compensation

When considering a wider area, i.e. decompacting a 2D transect or a 3D basin, a flexural compensation is required because a surface load tends to influence its surroundings. An Airy local compensation would not capture this effect and thus would overestimate the isostatic correction to be applied during the backstripping or backtracking (Roberts et al., 1998). Flexural compensation is based on the flexural strength of the lithosphere that is defined by its flexural rigidity D (in N m):

$$\begin{array}{}\text{(12)}& {\displaystyle }{\displaystyle }D={\displaystyle \frac{E{T}_{\text{e}}^{\mathrm{3}}}{\mathrm{12}(\mathrm{1}-{\mathit{\nu }}^{\mathrm{2}})}},\end{array}$$

where E corresponds to the Young modulus (N m⁻²), T_e is the effective elastic thickness of the lithosphere (m), and ν is the Poisson ratio. Total 1D flexure of the lithosphere for a line load on an infinite plate is given by the general analytical solution (Turcotte and Schubert, 2002):

$$\begin{array}{}\text{(13)}& {\displaystyle }{\displaystyle }D{\displaystyle \frac{{\mathrm{d}}^{\mathrm{4}}w}{\mathrm{d}{x}^{\mathrm{4}}}}+({\mathit{\rho }}_{\mathrm{m}}-{\mathit{\rho }}_{\mathrm{w}})gw=q,\end{array}$$

and its 2D expression

$$\begin{array}{ll}{\displaystyle }D{\mathrm{\nabla }}^{\mathrm{4}}w\left(x\right)& {\displaystyle }=D{\displaystyle \frac{{\mathrm{d}}^{\mathrm{4}}w}{\mathrm{d}{x}^{\mathrm{4}}}}+D{\displaystyle \frac{{\mathrm{d}}^{\mathrm{4}}w}{\mathrm{d}{y}^{\mathrm{4}}}}+\mathrm{2}-D{\displaystyle \frac{{\mathrm{d}}^{\mathrm{4}}w}{\mathrm{d}{x}^{\mathrm{2}}\mathrm{d}{y}^{\mathrm{2}}}}\\ {\displaystyle }& {\displaystyle }+({\mathit{\rho }}_{\mathrm{m}}-{\mathit{\rho }}_{\mathrm{w}})gw\left(x\right)\\ \text{(14)}& {\displaystyle }& {\displaystyle }=q,\end{array}$$

where w is the deflection of the lithosphere (m), g is the gravity acceleration (m s⁻²), q(x,y) is the vertical load (N m⁻²) applied to the lithosphere, and x and y are the coordinates in the horizontal plane. In this equation, downward-deflected lithosphere is substituted with seawater ρ_w. The load q(x,y) (in N m⁻²) associated with each ith sediment layer is calculated as follows:

$$\begin{array}{}\text{(15)}& {\displaystyle }{\displaystyle }q\left(i\right)={\mathit{\rho }}_{\text{ds}}\left(i\right)g{S}_{\mathrm{d}}\left(i\right).\end{array}$$

In PALEOSTRIP, 2D and 3D finite difference versions of Eqs. (13) and (14) have been implemented (e.g. Eqs. 9 and 10 in Wickert, 2016). They are based on Chapman (2015) flex2d and Cardozo (2009) flex3dv MATLAB codes that have been adapted to PALEOSTRIP needs. flex3dv routine is available at http://www.ux.uis.no/~nestor/Public/flex3dv.zip (last access: August 2021), and flex2d is available at https://www.jaychapman.org/matlab-programs.html (last access: August 2021). By means of a finite difference scheme, flexure can be computed by mean of an analytical solution (e.g. gflex, Wickert, 2016), whereas other numerical methods use the superposition of local solutions to point loads in the wavenumber domain (e.g. Green's functions, TAFI v1.0, Jha et al., 2017) or in the spectral domain (Wienecke et al., 2007), for example. These use biharmonic equation for plate flexure with uniform elastic properties. Approaches using the convolution method also allow the use of spatially variable elastic properties (e.g. Braitenberg et al., 2003, 2002). A file of spatially variable T_e can be provided through the GUI (Fig. 1a, bottom). A spatially uniform T_e can also be prescribed directly from the GUI. All other parameters involved in the flexural or Airy isostasy can be modified from the GUI (density constants, Poisson ratio, and Young modulus).

3.3 Thermal subsidence

Thermal subsidence corresponds to a vertical contraction of the lithosphere. During a rifting phase, as a first step, the lithosphere stretches apart and thins, which causes a net increase in heat outflow towards the surface due to the upwelling of the underlying asthenosphere. The stretching is generally not uniform, but reconstructions of past stretching processes require knowledge about plate tectonic strain rates (e.g. Müller et al., 2019), which is beyond the scope of our software. The first step is called “initial” subsidence. In the second step, the “thermal subsidence” occurs; i.e. the stretched lithosphere constricts due to cooling. To account for those effects during backstripping, PALEOSTRIP adopts the 1D-thermal subsidence model from McKenzie (1978) that assumes an instantaneous stretching (syn-rift) and a single rifting phase. The initial subsidence Z_init is given by

$$\begin{array}{}\text{(16)}& {\displaystyle }{\displaystyle }{Z}_{\text{init}}={\displaystyle \frac{T_{\text{IL}}\left[({\mathit{\rho }}_{\mathrm{m}}-{\mathit{\rho }}_{\mathrm{c}})\frac{T_{\text{IC}}}{T_{\text{IL}}}\left(\mathrm{1}-\frac{{\mathit{\alpha }}_v{T}_{\mathrm{m}}{T}_{\text{IC}}}{\mathrm{2}{T}_{\text{IL}}}\right)-\frac{{\mathit{\rho }}_{\mathrm{m}}{\mathit{\alpha }}_v{T}_{\mathrm{m}}}{\mathrm{2}}\right]\left(\mathrm{1}-\frac{\mathrm{1}}{\mathit{\beta }}\right)}{{\mathit{\rho }}_{\mathrm{m}}(\mathrm{1}-{\mathit{\alpha }}_v{T}_a)-{\mathit{\rho }}_{\mathrm{w}}}},\end{array}$$

where T_IL and T_IC are the initial lithospheric and crustal thicknesses at the beginning of rifting, β is the stretching factor, α_v is the coefficient of thermal expansion, ρ_c is the crustal density, and ρ_m is the mantle density. In this equation, the subsidence is isostatically compensated by water (ρ_w) using the Airy local compensation. This 1D instantaneous model can be improved by adopting a time-evolving approach to the extension, for example that of Jarvis and McKenzie (1980). However, comparison between Jarvis and McKenzie (1980) and McKenzie (1978) revealed that the two models show no or only little discrepancy if the duration of extension, given the time required to extend it by a factor of β, is about $\mathrm{60}/{\mathit{\beta }}^{\mathrm{2}}$ Myr if β≤2 or $\mathrm{60}(\mathrm{1}-\mathrm{1}/\mathit{\beta }{)}^{\mathrm{2}}$ Myr if β≥2 (Jarvis and McKenzie, 1980). During the second step, the thermal subsidence occurs after the end of rifting (post-rift) and accounts for the vertical thermal conduction, which takes the shape of an exponential function decaying with time:

$$\begin{array}{}\text{(17)}& {\displaystyle }{\displaystyle }{Z}_{\text{thermal}}\left(t\right)=E_{\mathrm{0}}{\displaystyle \frac{\mathit{\beta }}{\mathit{\pi }}}\sin \left({\displaystyle \frac{\mathit{\pi }}{\mathit{\beta }}}\right)\left(\mathrm{1}-\exp \left(-{\displaystyle \frac{t}{\mathit{\tau }}}\right)\right),\end{array}$$

where t corresponds to the time elapsed since the end of rifting expressed in seconds (e.g. time of backtracked horizon = 14 Ma; age of the end of rifting = 76 Ma; and $t=(\mathrm{76}-\mathrm{14})\times \mathrm{3600}\times \mathrm{24}\times \mathrm{365}\times {\mathrm{10}}^{\mathrm{6}}$).

$$\begin{array}{}\text{(18)}& {\displaystyle }{\displaystyle }{E}_{\mathrm{0}}={\displaystyle \frac{\mathrm{4}{\mathit{\rho }}_{\mathrm{m}}{\mathit{\alpha }}_v{T}_{\mathrm{m}}{T}_{\text{IL}}}{{\mathit{\pi }}^{\mathrm{2}}({\mathit{\rho }}_{\mathrm{m}}-{\mathit{\rho }}_{\mathrm{w}})}},\text{ and }\mathit{\tau }={\displaystyle \frac{T_{\text{IL}}^{\mathrm{2}}}{{\mathit{\pi }}^{\mathrm{2}}\mathit{\kappa }}},\end{array}$$

where κ is the thermal diffusivity. The total thermal subsidence is given by

$$\begin{array}{}\text{(19)}& {\displaystyle }{\displaystyle }{Z}_{\text{tot\_thermal}}=Z_{\text{init}}+Z_{\text{thermal}}.\end{array}$$

In PALEOSTRIP, during the backtracking procedure, paleo-water depths are corrected by removing the increment of thermal subsidence between each time step and present-day. This is because backtracking is an iterative process departing from a known state, i.e. present day. Because Z_init is constant in time it cancels out, and the thermal subsidence correction is given by

$$\begin{array}{}\text{(20)}& {\displaystyle }\mathrm{\Delta }{Z}_{\text{tot\_thermal}}\left(t\right)={\displaystyle }(Z_{\text{init}}+Z_{\text{thermal}}(t\left)\right)-(Z_{\text{init}}+Z_{\text{thermal}}(\mathrm{0}\left)\right),\end{array}$$

and thus,

$$\begin{array}{}\text{(21)}& {\displaystyle }{\displaystyle }\mathrm{\Delta }{Z}_{\text{tot\_thermal}}\left(t\right)=Z_{\text{thermal}}\left(t\right)-Z_{\text{thermal}}\left(\mathrm{0}\right),\end{array}$$

where t varies from present to past in this equation following backstripping procedure, i.e. t is an age older than present day (0). ΔZ_{tot_thermal} is thus positive, with Z_thermal being larger for younger times, i.e. a larger time is elapsed since the end of rifting.

Note that by applying the 1D model from McKenzie (1978) to 2D transects and 3D maps, it is assumed that no horizontal heat advection occurs. Almost all parameters involved in the thermal subsidence can be modified (age of end of rifting, T_IL and T_IC, α_v, κ) from the GUI tab interface (Fig. 1b). The stretching factor β can be formulated through several methods.

• (1)

  By prescribing a constant and uniform β factor that will be used at each time step of backtracking to calculate the thermal subsidence.

• (2)

  By linearly interpolating (in the x direction for 3D maps) between two prescribed constant and uniform stretching factors, β₁ and β₂

• (3)

  By inputting a user-based constant but spatially variable β factor in one direction (x grid) to compute spatially variable thermal subsidence for a 2D transect or by inputting a 2D (x,y grid) file to compute spatially variable thermal subsidence for 3D maps.

3.4 Sea Level correction

On long timescales, sea level has been varying with time due to plate tectonics changing the dimensions of ocean basins (e.g. Müller et al., 2018a) and, on shorter timescales, due to continental ice storage within ice sheets, ice caps, and glaciers during cold periods, as was the case, e.g. during the second half of the Cenozoic (the last 34 Myr, e.g. Miller et al., 2020). Classically, only eustatic sea level changes have been considered for correcting the subsidence or the paleo-water depth. For the last 34 Myr, eustatic sea level is usually defined relative to the present-day total ocean area (362.15×10⁶ km²) because it is assumed that this area evolved only a little and not enough to significantly alter this number. For time periods older than that, eustatic sea level changes have to be expressed according to the ocean area of the time. Considering only eustatic sea level changes results in highly approximated ice sheet growth and decay, leading to sea level variations of growing amplitude over the past 34 Myr (e.g. Miller et al., 2020). The variations induced by glaciations have much shorter timescales, i.e. 10–10⁵ years, compared with those induced by plate tectonics (10⁵–10⁶ years). Associated sea level changes are not spatially uniform and induce changes in the Earth's gravity field (e.g. Tamisiea et al., 2001), and regional self-gravitating sea level changes can substantially vary from eustasy (e.g. Tamisiea et al., 2001; Clark et al., 2002; Milne and Mitrovica, 2008).

In PALEOSTRIP, sea level (SL, in m) is corrected the same way as thermal subsidence:

$$\begin{array}{}\text{(22)}& {\displaystyle }{\displaystyle }\mathrm{\Delta }\text{SL}\left(t\right)=\text{SL}\left(t\right)-\text{SL}\left(\mathrm{0}\right),\end{array}$$

where t varies from present to past in this equation following a backtracking procedure, i.e., t is an age older than present day (0). ΔSL is expressed relative to present, and it can therefore assume positive or negative values, i.e. induce an uplift or a subsidence. Most of the sea level variation time series found in the literature already correspond to variations relative to present, i.e. are already in the form of Eq. (22). Thus, in order to avoid confusion for the user, three different ways of correcting water depth with sea level changes have been implemented within PALEOSTRIP and can be managed through the GUI (Fig. 1c):

• (1)

  by prescribing constant and uniform sea level correction ΔSL applied at each time step of the backtracking (e.g. correction at a time where sea level is lower by 100 m relative to present is prescribed as −100 in the GUI),

• (2)

  by applying a spatially uniform but time-varying sea level correction based on time series from Haq et al. (1987) or from Miller et al. (2020) (implemented within PALEOSTRIP);

• (3)

  by applying a user-provided time series or a constant in time but spatially varying map $(x,y,z)$ of sea level changes relative to present.

The last option allows us to prescribe sea level changes calculated with glacio-isostatic adjustment models (e.g. SELEN⁴ Spada and Melini, 2019) to account for regional self-gravitating effects of ice sheet growth and decay. ΔSL is used to compute the water load correction to adjust the water depth WD. In Eq. (3), water load due to sea level change is described based on Airy local compensation. In PALEOSTRIP, the water load due to sea level change is computed using the isostatic method previously selected to carry out decompaction (see Sect. 3.2). Note that if isostasy is deactivated, no water load correction is computed.

Figure 5 Close-up of the Ross Sea bathymetry from IBCSO (Arndt et al., 2013) and its location in Antarctica. Most of the ice-free bathymetry is backtracked as shown in the case studies below. The BGR80-007 2D marine seismic transect used for validation is indicated with a red line, and the Deep Sea Drilling Project (DSDP) 273 site location, backtracked in the case studies below, is indicated with a yellow square.

Figure 6 The 2D transect BGR80-007 case study. (a) Input present-day compacted depths (m) of bathymetry (yellow) and the nine seismic unconformities, including the mid-Miocene Ross Sea Unconformity 4 (RSU4) and basement (thick blue line). Layout is from the PALEOSTRIP Plot GUI. Comparison of backtracked basement depths for different but spatially uniform lithospheric elastic thicknesses (Te): (b) between the PALEOSTRIP finite difference isostasy and the Flex-Decomp (Kusznir et al., 1995) fast Fourier transform (wavenumber-based) isostasy model; (c) between PALEOSTRIP finite difference isostasy and the TAFI (Jha et al., 2017) Green's function spectral isostasy model. Panel (d) is the same as (b) but accounts for the thermal subsidence correction using a rift age of 85 Myr and a uniform β of 2.

Figure 7 PALEOSTRIP Application GUI plot and save results interface.

3.5 Dynamic topography

Mantle dynamics generate flows that cause time-varying surface topography and bathymetry deformations, which are called dynamic topography. The timescale at which the mantle flow produces dynamic topography that occurs at long wavelengths (≈5000 to 10 000 km Hoggard et al., 2016). The current mantle-driven dynamic topography was revealed by estimating the residual topography, obtained after removing the isostatic topography, which is generated by thickness and density contrasts within the lithosphere and the isostatic effect of the lithosphere calculated using seismic data compilations (e.g. Flament et al., 2013; Hoggard et al., 2016). Stratigraphic observations of continental margins revealed that the dynamic topography evolves quickly with time, potentially impacting long-term climate evolution (Hoggard et al., 2016). Austermann et al. (2017) showed that the magnitude of past interglacial sea level proxies partly results from dynamic topography, and they suggested correcting sea level proxies before inferring the relative contribution of past ice sheet to the sea level changes at that time. Furthermore, simulated Pliocene dynamic topography changes accounted for in Antarctic ice sheet simulations revealed that ice sheet stability is highly influenced by mantle dynamics that create or cancel pinning areas at the surface (Austermann et al., 2015). Thus past reconstructions of continental morphology and shallow margins should account for past evolution of dynamic topography. Müller et al. (2018a) recently provided time slices of reconstructed dynamic topography over the past 240 Myr, which constitutes a strong basis to calculate a dynamic topography correction.

PALEOSTRIP provides three ways to correct for dynamic topography changes through its GUI (Fig. 1d):

1. by prescribing a uniform and constant dynamic topography change relative to present day that will be applied to each time step of the backtracking procedure;

2. by using a user-based spatially uniform time series of dynamic topography changes relative to present;

3. by using user-based 2D maps of dynamic topography changes relative to present, e.g. Müller et al. (2018a) (note that Müller et al., 2018a, is not implemented within PALEOSTRIP); maps of dynamic topography are inputs to PALEOSTRIP and the user is free to use any reconstructions (note that inputs of dynamic topography require some pre-processing to be adjusted to the area of interest before being passed through the GUI).

The correction is calculated as for sea level changes (Eq. 22) as follows:

$$\begin{array}{}\text{(23)}& {\displaystyle }{\displaystyle }\mathrm{\Delta }\text{DynT}\left(t\right)=\text{DynT}\left(t\right)-\text{DynT}\left(\mathrm{0}\right),\end{array}$$

where t varies from present to past in this equation following backstripping procedure, i.e., t is an age older than present (0). ΔDynT is expressed relative to present, and since dynamic topography is not spatially uniform at a global level, it can therefore assume positive or negative values, i.e. induce an uplift or a subsidence.

3.6 Sediment erosion

Sediment erosion and depositions are the main processes involved in the building of continental margins. Present-day identified sedimentary units may not reflect the real number of sedimentary units at a given time in the past because some of the layers might have been eroded in the meantime. Erosion is identified within sediment cores as hiatus and as high-amplitude (sometimes truncational) reflectors within marine seismic profiles. During backstripping, one should move back eroded sediments to their presumed original locations (e.g. Paxman et al., 2019), which can be either within the backstripped area or outside the backstripped area. If eroded sediments are moved back in the backstripped area, then the total number of layers should be modified during the backward modelling process. The current version of PALEOSTRIP does not treat erosion and does not allow for the variation of the total number of sedimentary layers once they are input at the beginning of the procedure. A good description of how to handle erosion during backstripping is provided in Sect. 5.4 of Wangen (2010).

4 PALEOSTRIP grid interpolation

PALEOSTRIP handles 1D data (drill sites), 2D data (transects), and 3D data (grids). All calculations (except flexure) are performed on the grid or array of original input data. For 2D and 3D data, all input horizon depths have to be on the same grid or array (identical coordinates), and spacing along the X and Y directions must be constant. For 3D grids, spacing in the X direction can differ from spacing in the Y direction. For example, 2D transects must have horizons depths of the same length. The 3D maps can be provided as an irregular polygon, and all horizon depths must be provided on the same exact polygon (see examples in Sect. 6). At last, grid spacing dx and dy must be integers (e.g. dx=10 km and not dx=4.3 km).

4.1 2D data: vertical transects

Vertical transects imply that input data are provided along an horizontal direction X and a vertical direction (depth) Z. Due to the needs of flexure calculations, the original sediment loads array is extended (duplication of last values at both edges) to about 90 % of its length from both edges (Fig. 4a). The sediment loads are then placed at the centre of an expanded domain (3 times larger than the original array length) to avoid edge effects on flexure correction. After flexure calculation, the correction is extracted and relocated on the original array domain. All the other backtracking calculations occur on the original input array. Note that spatially variable lithospheric elastic thickness is also interpolated and extrapolated following this procedure.

4.2 3D data: maps

Maps imply that data are provided along the two horizontal directions, X and Y, and along the vertical direction, Z. Input data are read as scattered data and not as gridded data. This means that points are unstructured and are not written in a file following a classical NX×NY structure but are instead treated as independent single points. When data are treated as scattered this takes more computational time, but this also allows one to input either structured gridded data or irregular polygon data. This facilitates all computations and avoids unnecessary duplication of routines for 1D, 2D, or 3D cases. For the need of flexure, 3D data are interpolated (preserving their original horizontal resolution) on a regular rectangular grid. The original input grid is expanded (extrapolation of the last values from the edges) to about 30 % from all sides of the grid. The sediment loads are then placed at the centre of an expanded domain (twice as large as the original grid dimension) to avoid edge effects on flexure correction (Fig. 4b). After flexure calculation, the correction is extracted and relocated on the original grid domain. All the other backtracking calculations occur on the original input grid. Note that spatially variable lithospheric elastic thickness is also interpolated and extrapolated following this procedure.

Figure 8 Backtracking steps for the western Ross Sea DSDP 273 well site. Present-day depths of the basement (black), Ross Sea unconformity 5 (blue), Ross Sea Unconformity 4A (orange), and bathymetry (brown) at the time of RSU4A (19 Ma), RSU5 (21 Ma), and at the time before glacial sediment deposition (>34 Ma). For each time step, two backtracking are shown, one accounting for Airy isostatic correction only and one accounting for Airy isostatic correction and thermal subsidence (dashed). Physical parameter values used are displayed in Table A1 in Appendix A.

Figure 9 Initial input present-day depths (m): (a) bathymetry, (b) mid-Miocene unconformity RSU4, and (c) basement both from ANTOSTRAT atlas (Brancolini et al., 1995). Layout is from the PALEOSTRIP plotting GUI; the colour scale is the default colour scale implemented within PALEOSTRIP. The colour scale has been saturated below −4000 m and above 200 m for this figure: the maximum depth of the basement reaches −10 km and would have largely damped the other bathymetric changes occurring at shallower depths in panels (a, b). The island located in the uppermost right corner is the only location with an elevation above sea level. Note that the colour scale is colour-blind friendly and is from Crameri et al. (2020). PALEOSTRIP has implemented 12 different colour scales, and 9 out of 12 are colour-blind friendly.

Figure 10 Backtracked mid-Miocene bathymetries at about 14 Ma: (a) corrected only for Airy isostasy, (b) the same as (a) but including thermal subsidence correction with a spatially uniform β=2.1 and the end of rifting sets at 76 Ma, (c) the same as (b) but including eustatic sea level correction based on Miller et al. (2020) δ¹⁸O-derived reconstruction, (d) the same as (c) but including dynamic topography correction based on the Müller et al. (2018a) geodynamical model M1 (following Hochmuth et al., 2020). Layout is from the PALEOSTRIP plotting GUI; the colour scale is the default colour scale implemented within PALEOSTRIP and has been saturated below −4000 m and above 200 m for this figure. The island located in the uppermost right corner is the only location with an elevation above sea level. Note that the colour scale is colour-blind friendly and is from Crameri et al. (2020). PALEOSTRIP has implemented 12 different colour scales, and 9 out of 12 are colour-blind friendly. For all parameters used in these examples, see Table A1.

Figure 11 Backtracked mid-Miocene bathymetries at about 14 Ma. The same as Fig. 10 but using flexural isostasy with spatially variable lithospheric effective elastic thickness from Chen et al. (2018). Layout is from the PALEOSTRIP plotting GUI; the colour scale is the default colour scale implemented within PALEOSTRIP and has been saturated below −4000 m and above 200 m for this figure. The island located in the uppermost right corner is the only location with elevation above sea level. Note that the colour scale is colour-blind friendly and is from Crameri et al. (2020). PALEOSTRIP has implemented 12 different colour scales, and 9 out of 12 are colour-blind friendly. For all parameters used in these examples, see Table A1 in Appendix A.

5 PALEOSTRIP validation

We backtrack a 2D transect with PALEOSTRIP and with Flex-Decomp (Kusznir et al., 1995) to validate the results. The transect used in the case study is a revised version of the one studied by De Santis et al. (1999): the BGR80-007 seismic profile. The transect BGR80-007 is composed of nine identified seismic stratigraphic unconformities. The transect is about 250 km long, is broadly oriented north–south, and is located in the eastern Ross Sea (Fig. 5). The set of initial data is composed of 10 files: the actual bathymetry of the Ross Sea and the present-day depth of the nine seismic unconformities, including the present-day depth of the basement below (Fig. 6a). Depths are related to present regional sea level.

The 11th file contains the lithological parameters of the layers to be decompacted, as well as other parameters needed by PALEOSTRIP, excluding present-day bathymetry: LAYER is the layer number (1 to N), from bottom to top), POROSITY is the deposition porosity (unitless), DEC CON (1 km⁻¹) corresponds to the porosity decompaction coefficient, MAT DEN (kg m⁻³) corresponds to the compacted sediment density, AGE BASE (millions of years ago, Ma) is the age of horizons at the base of the layers, and NAME is the string name of horizons. The parameters are taken from the study of De Santis et al. (1999) and the lithological parameter file uses the following format:

NUMBER OF LAYERS =           9

LAYER POROSITY   DEC CON  MAT DEN    AGE BASE NAME
                 (1/KM)   (KG / M3)   (MA)

  1    0.4900    0.2700    2680.00    95.00   basement
  2    0.4500    0.4500    2680.00    26.00   rsu_6
  3    0.4500    0.4500    2680.00    24.90   rsu_5b
  4    0.4500    0.4500    2680.00    19.70   rsu_5a
  5    0.4500    0.4500    2680.00    18.00   rsu_5
  6    0.4500    0.4500    2680.00    14.20   rsu_4
  7    0.4500    0.4500    2680.00    10.00   rsu_3
  8    0.4500    0.4500    2680.00    4.00    rsu_2
  9    0.4500    0.4500    2680.00    0.60    rsu_1

To facilitate the comparison, thermal subsidence, sea level, and dynamic topography are switched off both in PALEOSTRIP and in Flex-Dedcomp. We compare the final backtracked depths of the basement using Airy local isostasy and flexural isostasy with different lithospheric elastic thicknesses (Fig. 6b). The match between PALEOSTRIP and Flex-Decomp is very good and persistent discrepancies (a few tens to hundreds of metres) are likely due to (1) the different ways of computing flexural isostasy in the spectral domain for Flex-Decomp and with finite difference for PALEOSTRIP, (2) re-interpolation of the load to a different resolution in Flex-Decomp (no reinterpolation in PALEOSTRIP), and (3) different extrapolation of the load outside of the original transect length to avoid edge effects on flexure. In PALEOSTRIP, the last point of the transects at both edges is duplicated to extend the original length of about 80 % at both sides (Fig. 4a). We also compare PALEOSTRIP backtracked results with those using the analytic solution from TAFI v1.0 (Jha et al., 2017) and implemented within PALEOSTRIP for the need of comparison (Fig. 6c). Similarly to the comparison with Flex-Decomp, we only account for isostasy and the other processes are switched off. Re-interpolation of the loads is performed by PALEOSTRIP in both cases. Comparison between PALEOSTRIP (finite difference scheme) and TAFI v1.0 reveal almost identical results.

Finally, we test the isostasy and thermal subsidence model by comparing those obtained by PALEOSTRIP and Flex-Decomp. De Santis et al. (1999) originally used Flex-Decomp with β=2, with the age of rifting set to 85 Ma and various lithospheric elastic thickness values to restore paleo-bathymetries. We use the same parameters on this revised transect in both Flex-Decomp and PALEOSTRIP. Match between PALEOSTRIP and Flex-Decomp is very good (Fig. 6d).

6 Case study: example of the Ross Sea

PALEOSTRIP GUI presents a plot and save interface to support each step of backtracking: the user can plot initial input data, backtracked data, and calculated intermediate variables relevant to the backtracking process and save them to ASCII files (Fig. 7). The user can also extract some 2D transects or 1D wells from 3D backtracked maps or 2D transects and plot and save them to ASCII files. In the following, we provide two case studies to illustrate the possibilities of PALEOSTRIP. Note that the physical parameters, sea level correction, or thermal-subsidence-related variables are not tuned since the aim of the examples is to illustrate the physics of PALEOSTRIP rather than to provide a realistic reconstruction of the paleo-bathymetries of this area.

Both the cases are taken from the continental margin in the western Pacific sector of Antarctica in the Ross Sea (Fig. 5). Sediment layers mostly accumulated after the main rifting phases that occurred in this area between 95 and 79 Ma (Stock and Cande, 2002; Decesari et al., 2007; Siddoway et al., 2004). Ross Sea stratigraphic data are currently being revised and differ from the ANTOSTRAT data (Brancolini et al., 1995) upon which the following 3D example is based. New reconstructed Ross Sea paleo-bathymetry using revised data will be the object of a specific contribution. ANTOSTRAT data have been recently used in pan-Antarctic reconstructions of past topographies and bathymetries (Paxman et al., 2019; Hochmuth et al., 2020).

6.1 Well (1D): Deep Sea Drilling Project (DSDP) site 273 – Ross Sea

In this example, we decompact the drilling site DSDP273 Hayes and Frakes (1975) from the western Ross Sea (Fig. 5). It has two main identified seismic unconformities above the basement. Lithological parameters are taken from De Santis et al. (1999), and the lithological parameters file (see Sect. 5) is as follows:

NUMBER OF LAYERS =           3

LAYER POROSITY  DEC CON   MAT DEN     AGE BASE NAME
                 (1/KM)   (KG / MC)   (MA)

  1    0.4500    0.2700    2680.00    95.00   basement
  2    0.4500    0.4500    2680.00    21.00   rsu_5
  3    0.4500    0.4500    2680.00    19.00   rsu_4a

Backtracking is performed using Airy local isostasy since flexure cannot be applied to 1D drilling sites. We also add thermal subsidence in order to illustrate the difference in backtracked depths of those unconformities (Fig. 8). The layout of Fig. 8 is not the layout of PALEOSTRIP, but the results have been assembled to highlight the impact of thermal subsidence on backtracked depths. As for 2D transects, all intermediate variables and input conditions can be plotted and saved.

6.2 Map (3D): ANTOSTRAT data

In this example we backtrack the depth of the mid-Miocene unconformity (14 Ma) across the Ross Sea from ANTOSTRAT atlas (Brancolini et al., 1995). The initial grid is an irregular polygon. The set of initial data is composed of three files: the actual bathymetry of the Ross Sea, the present-day depth of mid-Miocene unconformity, and the presumed present-day depth of the basement below (Fig. 9).

The fifth file contains the lithological parameters of the layers to be decompacted (see Sect. 5) and are taken from De Santis et al. (1999). The lithological parameters file is as follows:

NUMBER OF LAYERS =           2

LAYER  POROSITY  DEC CON   MAT DEN    AGE BASE  NAME
          (%)    (1/KM)   (KG / CC)   (MA)

  1    0.450     0.450     2680.0     95.000    Basement
  2    0.450     0.450     2680.0     14.200      RSU4

PALEOSTRIP is run several times to add one of the following components at a time: isostasy, thermal subsidence, sea level correction, and dynamic topography correction. Thanks to this approach, the user can perform ensembles of simulations to retrieve sound statistics of the model parameter space. Two series are shown, one using the Airy isostatic correction (Fig. 10) and the other using the flexural isostatic correction (Fig. 11). Paleo-bathymetry retrieved with flexural isostasy produces a quite different morphology from the one computed using Airy isostasy, as already observed by Roberts et al. (1998). Sea level and dynamic topography (Fig. A1b) corrections are not big enough to produce a significant change of the overall morphology. However, they matter for the bathymetric highs (especially the shallowest one), as a few tens of metres can uplift those highs above sea level. PALEOSTRIP allows the user to plot different variables amongst initial input data and computed quantities, such as isopach, density, porosity, or isostatic correction, allowing the user to check and separate the various processes required to perform a detailed analysis of their impact on the paleo-bathymetric reconstruction (Fig. A1).

7 Conclusions

PALEOSTRIP is one of the first examples of open-source 3D backtracking software. It can process paleo-bathymetries for 1D drilling sites, 2D transects, and 3D maps. It allows users to separate the various processes involved in the backtracking procedure. Thanks to this approach, the user can perform ensembles of simulations to retrieve sound statistics about the model parameter space. PALEOSTRIP has been designed to be modular to allow users to insert their own modifications. The code is documented, and thus implementation of new modules would only require minor work.

Appendix A: Case studies: settings

A1 Table of physical parameters

In all the examples provided, we used PALEOSTRIP default values automatically inserted within GUI (Table A1).

Miller et al. (2020)Müller et al. (2018b)

Table A1 Physical parameters values used in the case studies.

A2 Backtracked intermediate variables

Here we plot some of the intermediate physical variables calculated during the backtracking procedure and available for plotting through PALEOSTRIP GUI: sediment isopachs, decompacted density, decompacted porosity, isostatic correction, and dynamic topography (Fig. A1). Thermal subsidence is also available for plotting, but since we employ a spatially uniform β value here, thermal subsidence is consequently spatially uniform over the domain, and thus we do not display it.

Figure A1 Example of input data and related backtracked sediment variables during the mid-Miocene: (a) spatially variable lithospheric elastic thickness (Chen et al., 2018) and (b) dynamic topography correction (Müller et al., 2018b). Both (a) and (b) are user-provided input to PALEOSTRIP. (c) Flexural isostatic correction calculated using (a). (d) Decompacted isopach (m). (e) Decompacted porosity (%). (f) Decompacted density (kg cm⁻³). Note that the colour scale is colour-blind friendly and is from Crameri et al. (2020). PALEOSTRIP has implemented 12 different colour scales, and 9 out of 12 are colour-blind friendly.