Tectonic forcing

Note

These forcing are user defined and could be either lithospheric or mantle induced but goSPL does not know which underlying process is inducing these changes and does not account for sediments and crust compression.

Uplift & subsidence

In the most simple case of vertical-only displacements (i.e. uplift or subsidence), the model requires the declaration of successive events using start, end times and the tectonic rates map upsub (in metres per year). This map is defined prior to the model run for goSPL spherical or 2D mesh as a numpy zip array.

During the model run and for each time step, the prescribed rates are applied on all vertices of all partition by increasing or decreasing their elevation values following imposed local subsidence or uplift rates.

Horizontal advection

The second, more complex, option consists in adding horizontal displacements (hdisp). In this case, the user has to defined for each point of the grid the associated velocities (represented by displacement fields in both X,Y and Z coordinates) and advection could be performed based on different techniques.

Using Upwind or IIOE scheme

Important

These two approaches are designed for 2D meshes and have not been tested on global simulations.

The Upwind scheme offers a rapid solution to the advection problem but with potentially excessive diffusion when solved implicitly (first-order upwind implicitly scheme).

Alternatively, the user can perform the advection of elevation and erosion by solving the advection equation based on a finite volume space discretization and a semi-implicit discretization in time called the Inflow-Implicit/Outflow-Explicit scheme and following the work from Mikula & Ohlberger, 2014.

Setting the approach for a given simulation is done in the input file by setting the advect key to either upwind, iioe1 or iioe2.

Note

In the IIOE approach, the outflow from a cell is treated explicitly while inflow is treated implicitly.

Since the matrix of the system is determined by the inflow fluxes it is an M-matrix yielding favourable solvability and stability properties.

Both methods allow large time steps without losing stability and not deteriorating precision.

The IIOE is formally second order accurate in space and time for 1D advection problems with variable velocity and numerical experiments indicate its second order accuracy for smooth solutions in general.

Note

Velocity at the face is taken to be the linear interpolation for each vertex (in a vertex-centered discretisation the dual of the delaunay triangulation (i.e. the voronoi mesh has its edges on the middle of the nodes edges)).

Similarly we consider that the advected variable at the face is defined by linear interpolation from each connected vertex.

Two IOOE schemes are available the first one (iooe1) uses the technique defined here, the second (iooe2) is more computationally demanding but should limit instabilities in regions with large advection velocities (see paper).

Using the semi-lagrangian scheme

Used for global mesh, the horizontal advection (advect: interp) often corresponds to plate motions and can be obtained from plate-tectonic reconstruction software.

Warning

Compared to vertical displacements applied at each time step, the horizontal advection for this method is only performed at the end of the specified tectonic interval (so usually of the order of several thousands to millions years).

Fig. 6 Illustration of interpolation between mesh vertices induced by horizontal velocities based on a semi-lagrangian scheme.

Note

Due to tectonic advection, the density of the surface nodes evolves over time, which leads to areas showing rarefaction or accumulation of nodes. In order for the finite volume scheme to remain accurate, local addition and deletion of nodes and remeshing of the triangulated spherical surface are therefore required. However, the remeshing process is computationally expensive and requires to rebuild not only the delaunay mesh but also the associated voronoi one and to redistribute the mesh and its vertices parameters with PETSc.

To avoid remeshing, the initial mesh is assumed to remain fixed and the advected points are then used to interpolate the advected nodes informations (elevation, erosion, deposition, stratigraphic information…) onto the fixed ones using an inverse weighting distance approach and SciPy cKDTree. The number of points used for the interpolation is defined by the user in the input file (interp field). This is an efficient approach but it might not be the most accurate one.

Flexural isostasy

The flexural isostasy in goSPL is computed on a regular grid, and therefore required to perform a series of interpolation from and to the goSPL unstructured mesh.

Fig. 7 Flexural isostasy can be produced in response to a range of geological loads (from Wickert, 2016).

gFlex for 2D simulations

When running goSPL in 2D, it is possible to compute the flexural isostasy equilibrium based on topographic change. The function accounts for flexural isostatic rebound associated with erosional loading/unloading using gFlex.

Important

The approach is performed in serial and therefore can be relatively slow depending on the simulation size.

It takes an initial (at time \(t\)) and final topography (at time \(t + \Delta t\)) (i.e. before and after erosion/deposition) and returns a corrected final topography that includes the effect of erosional/depositional unloading/loading.

The approach solves the bi-harmonic equation governing the bending/flexure of a thin elastic plate floating on an inviscid fluid (the asthenosphere).

\[D \frac{d^4 w}{d^4 x} + \Delta \rho g w = q\]

where \(D\) is the flexural rigidity, \(w\) is vertical deflection of the plate, \(q\) is the applied surface load, and \(\Delta \rho = \rho_m − \rho_f\) is the density of the mantle minus the density of the infilling material.

pyshtools for global simulations

When running goSPL in global mode, the flexural isostasy equilibrium is computed based on topographic change using pyshtools.

In this case, the global flexural-isostatic response is computed via thin-elastic-shell theory in the spherical-harmonic domain (Turcotte 1979; Willemann & Turcotte 1981).

There are two modes available for the flexural isostasy computation:

1. In the case of a constant Te (scalar te), a single-pass spectral solve is used.

2. In the case of spatially-varying Te (te), Picard iteration with the constant-Te spectral solver as the inner kernel; the rigidity perturbation D(x) - D0 is treated as a spatial correction load that is updated each iteration until self-consistent.

Note that the approximation drops gradient-of-D terms in the elastic operator, so very sharp Te steps may need to be smoothed (e.g. Gaussian with sigma ~ flexural wavelength / 4) before being passed in.