Tectonic related parameters#

Tectonic forcing parameters#

Declaration example:

tectonics:
    - start: -20000000.
      end: -19000000.
      upsub: ['data/uplift20Ma','t']
    - start: -18000000.
      end: -17000000.
      upsub: ['data/uplift18Ma','t']
      hdisp: ['data/hdisp18Ma', 'hxyz']

It defines the tectonic forcing conditions from a sequence of events defined by a starting and ending time (start and end) and either a vertical rate only forcing (e.g. uplift and/or subsidence defined with upsub) or a fully 3D displacement rate hdisp. These displacement rates are set in metres per year.

Important

For horizontal advection the user needs to define the advect key in the domain section of the input file. The advection scheme could either be upwind, iioe1, iioe2 or interp (go to the technical information in the documentation for more information).

Important

Here again, these forcing files are defined as numpy zip array (.npz). These files use specific keys to identify the tectonic forcing that are specified by the user in the input file. For vertical only condition, the displacements (in m/yr) is a 1D vector with values on each node of the grid. For the horizontal condition, the key is a 3D array containing the displacement rates along the x, y and z axis (in m/yr). When the horizontal advection is for 2D grids, the provided displacements also need to be in 3D but with the third dimension set to 0.0.

Note

There is no requirement to impose continuous tectonics forcing and you might chose to have periods without displacement by making discontinuous events using the start and end keys.

Note

When applying horizontal displacement using the interp scheme (mainly used in global simulation to represent plate movements), the horizontal movements are performed at the end of the period (i.e., at the specified end time). In the other cases, the horizontal displacement is done at every timestep (specified by dt).

Flexural isostasy definition#

This function computes the flexural isostasy equilibrium based on topographic change. It accounts for flexural isostatic rebound associated with erosional loading/unloading, solving the thin-elastic-plate biharmonic equation \(D\,\nabla^4 w + \Delta\rho\,g\,w = q\) with \(D = E\,T_e^3/[12(1-\nu^2)]\).

goSPL provides two solvers, selected by method:

'fem' (flat / 2D models, default) — a parallel finite-volume biharmonic solve directly on the unstructured mesh: no regular grid, no gather-to-root, no external dependency, and spatially-varying elastic thickness handled in a single linear solve. This is the only flat-model solver (it replaced the former 'FD' and 'FFT' approaches; those keys are no longer accepted).
'global' — a spherical-harmonic solve for global (sphere) models, run serially on the root rank.

Declaration example (flat / 2D model):

flexure:
  method: 'fem'
  thick: 30.e3
  rhoc: 2300.0
  rhoa: 3300.0
  young: 65e9
  nu: 0.25
  bcN: "Mirror"
  bcE: "0Slope0Shear"
  bcS: "Mirror"
  bcW: "0Slope0Shear"

where:

method: 'fem' for flat/2D models (default) or 'global' for global (spherical) models.
thick effective elastic plate thickness \(T_e\) in m (uniform; spatially/temporally variable \(T_e\) is set with the temap key below).
rhoc sediment/crust (load) density in kg/m³,
rhoa asthenospheric (mantle) density in kg/m³,
young Young’s Modulus in Pa,
nu Poisson ratio,
bcN, bcE, bcS, bcW North, East, South and West boundary conditions ('fem' only — see the boundary-condition note below).

The following keys apply to both methods:

interval: apply flexure only every interval goSPL steps, accumulating the load in between (default 1 = every step). The elastic response is instantaneous, so this only coarsens the temporal resolution of the loading history.

For the 'global' method (spherical-harmonic solve) the following additional keys are available:

i. ninterp: number of source points for the mesh ↔ Driscoll–Healy grid interpolation (default 4; unused by 'fem'). j. res_deg: resolution (degrees) of the Driscoll–Healy grid (default 0.25). The deflection is long-wavelength, so a coarser grid (e.g. 0.5 or 1.0) is markedly cheaper with little change to the result — the cheapest speed-up for the serial global solve. k. maxIter: maximum Anderson/Picard iterations for the varying-Te (temap) iterative solve (default 50). l. tol: relative convergence tolerance |dw|/|w| for that iteration (default 5.e-4). m. relax: under-relaxation factor (0–1) for the varying-Te update (default 1.0; lower damps strong Te contrasts).

Note

The varying-Te global solve warm-starts each step from the previous step’s converged deflection, so after the first step it typically needs noticeably fewer iterations. With a constant thick (no temap) the global solve is a single spectral pass and maxIter/tol/relax are unused.

Note

For flat/2D ('fem') models the per-side boundary conditions bcN/bcE/bcS/bcW (default 0Slope0Shear) can be:

0Slope0Shear: zero slope and zero shear at the edge (∂w/∂n = 0, ∂³w/∂n³ = 0) — the natural boundary; use it (or Mirror) when the load is far from the edge.

Mirror: reflective/symmetry boundary. For a thin plate this is identical to 0Slope0Shear and uses the same (natural) treatment.

0Displacement0Slope: clamped edge (w = 0 and ∂w/∂n = 0).

0Moment0Shear and Periodic are not available with the 'fem' solver.

Note

The 'fem' solver caches its operator and factorisation and reuses them every step (serial direct LU / parallel MUMPS), so it is fast even on small meshes; only the surface load changes between steps. Its result is influenced by the boundary condition (see bcN/bcE/bcS/bcW below) only where the flexural wavelength approaches the domain size (the deflection reaches the edges); away from the edges the interior solution is insensitive to the choice.

In addition, it is possible to define a spatially-variable (and time-varying) lithospheric elastic thickness with the temap key below, where specific maps can be added through time during the run. This works with both 'fem' and 'global' methods (for 'fem' the varying \(T_e\) is still a single linear solve).

Declaration example:

temap:
  - start: 250.e6
    map: ['input/temap1', 'te']
  - start: 251.e6
    map: ['input/temap2', 'te']

Here again, the elastic maps files are provided as numpy zip array (.npz).