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 imple-mentations and strengthens research synergies. To ﬁll 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 modiﬁed from the GUI. It includes 3D ﬂexural isostasy, 1D thermal subsidence, and possibilities to correct for prescribed sea level and dynamical topography following, we detail the physics embedded within PALEOSTRIP, and we show its application using a drilling site (1D), a transect 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.


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, 2013IPCC, , 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 2 , 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; Published by Copernicus Publications on behalf of the European Geosciences Union. 5286 F. Colleoni et al.: PALEOSTRIP 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(Lee et al., , 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. (2015Dressel et al. ( , 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.

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 MAT-LAB (2019) and runs on any operating system. Most of the code should be compatible with previous ver- sions 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.

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 PALE-OSTRIP 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.

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 . 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).

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. PALE-OSTRIP 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.

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): where the first term of the right-hand side corresponds to the isostatic compensation (ρ s sediment density, ρ m mantle den- sity, ρ 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 equa- . 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).
tion. 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: 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. socalled "dynamic topography" (Müller et al., 2018a). Paleowater depth needs to be corrected for dynamic topography changes DynT(i) causing uplift or subsidence of the topography and bathymetry: 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.

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: where φ 0 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 −1 , and z is depth in kilometres. This empirical equation implies that φ 0 is decreased by the factor of 1/e at the depth of 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: which upon integration gives where z i and z i+1 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: Based on the fact that V total = V sed + V 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: PALEOSTRIP solves this equation iteratively (starting with z i = 0 and then the further iterations) until convergence is reached. Lithological parameters, i.e. the decompaction coefficient c and the depositional surface porosity φ 0 , are prescribed in the user-provided parameter files (see Sect. 6).

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 and the decompacted sediment bulk density ρ ds (i) of the ith sediment layer is given by 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 PALE-OSTRIP, two methods are implemented.

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: 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).

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): where E corresponds to the Young modulus (N m −2 ), 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): and its 2D expression where w is the deflection of the lithosphere (m), g is the gravity acceleration (m s −2 ), q(x, y) is the vertical load (N m −2 ) 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 −2 ) associated with each ith sediment layer is calculated as follows: 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., 2003Braitenberg et al., , 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).

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 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 60/β 2 Myr if β ≤ 2 or 60(1 − 1/β) 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: 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 = (76 − 14) × 3600 × 24 × 365 × 10 6 ).
where κ is the thermal diffusivity. The total thermal subsidence is given by In PALEOSTRIP, during the backtracking procedure, paleo-water depths are corrected by removing the increment of thermal subsidence between each time step and presentday. 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 and thus, 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, β 1 and β 2 (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.

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 6 km 2 ) 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 5 years, compared with those induced by plate tectonics (10 5 -10 6 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: 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 4 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 PALE-OSTRIP, 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.

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: 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.

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).

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 sedi-ment 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.

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 N X × N Y 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 irregu- lar 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.

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 base-  Table A1 in Appendix A. ment 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 −1 ) corresponds to the porosity decompaction coefficient, MAT DEN (kg m −3 ) 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)  To facilitate the comparison, thermal subsidence, sea level, and dynamic topography are switched off both in PALE-OSTRIP and in Flex-Dedcomp. We compare the final backtracked depths of the basement using Airy local isostasy  (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.  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. 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 PALE-OSTRIP, (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 PA-LEOSTRIP, 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).  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. 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 PALE-OSTRIP. 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).

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). 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.

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 presentday 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: 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).

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.

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