Class GridProcess#

class tools.addprocess.GridProcess[source]#

When running goSPL on a 2D grid (i.e. not a global simulation), this class defines two processes operating on a regular grid:

Flexural isostasy: it allows to compute isostatic deflections of Earth’s lithosphere with uniform or non-uniform flexural rigidity. Evolving surface loads are defined from erosion/deposition values associated to modelled surface processes.
Orographic rain: it accounts for change in rainfall patterns associated to change in topography. The orographic precipitation function is based on Smith & Barstad (2004) linear model.

For global simulation, the library pyshtools provides a framework to estimate global-scale flexural isostasy.

Methods

`applyFlexure`()	This function computes the flexural isostasy equilibrium based on topographic change.
`cptOrography`()	Linear Theory of Orographic Precipitation following ` <https://journals.ametsoc.org/view/journals/atsc/61/12/1520-0469_2004_061_1377_altoop_2.0.co_2.xml>`_.

Initialise

__init__()

Initialisation of the gridProcess class.

Public Methods

`applyFlexure`()	This function computes the flexural isostasy equilibrium based on topographic change.
`cptOrography`()	Linear Theory of Orographic Precipitation following ` <https://journals.ametsoc.org/view/journals/atsc/61/12/1520-0469_2004_061_1377_altoop_2.0.co_2.xml>`_.

Private Methods

`_regInterp`(field)	Perform bilinear interpolation of `field` on the regular grid to unstructured 2D mesh.
`_buildRegGrid`()	Builds the regular grid based on nodes coordinates and instantiates two interpolation objects.
`_updateTe`()	Finds current elastic thickness map for the considered time interval.
`_cptFlex2D`(dZ)	Compute the flexural response for 2D cases.
`_buildDHGrid`()	Build the Driscoll-Healy (DH2, shape N x 2N) target grid for global flexure and precompute the interpolation weights linking the unstructured spherical mesh (`self.mCoords`) and the DH grid.
`_unstr2dh`(field)	Inverse-distance interpolation of a 1-D unstructured field defined at `self.mCoords` onto the DH2 grid (shape `(dh_N, 2*dh_N)`).
`_dh2unstr`(field_dh)	Bilinear interpolation of a DH2-grid array back to the unstructured mesh.
`_cmptFlexGlobal`(erodep_dh, te_dh[, ...])	Spectral thin-elastic-shell flexure solve on the precomputed Driscoll-Healy (DH2) grid.

Public functions#

GridProcess.applyFlexure()[source]#

This function computes the flexural isostasy equilibrium based on topographic change. It is a simple routine that accounts for flexural isostatic rebound associated with erosional loading/unloading.

The function takes an initial (at time t) and final topography (at time t+Dt) (i.e. before and after erosion/deposition) and returns a corrected final topography that includes the effect of erosional/depositional unloading/loading. It uses a spectral method to solve the bi-harmonic equation governing the bending/flexure of a thin elastic plate floating on an inviscid fluid (the asthenosphere).

\[D (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.

GridProcess.cptOrography()[source]#

Linear Theory of Orographic Precipitation following ` <https://journals.ametsoc.org/view/journals/atsc/61/12/1520-0469_2004_061_1377_altoop_2.0.co_2.xml>`_.

The model includes airflow dynamics, condensed water advection, and downslope evaporation. It consists of two vertically-integrated steady-state advection equations describing: (i) the cloud water density and (ii) the hydrometeor density. Solving these equations using Fourier transform techniques, derives a single formula relating terrain and precipitation.

Note

Please refer to the original manuscript of Smith and Barstad (2004) to understand the model physics and limitations.

Common bounds:

latitude : 0-90 [degrees]

precip_base : 0-10 [mm/h]

precip_min : 0.001 - 1 [mm/h]

conv_time : 200-2000 [s]

fall_time : 200-2000 [s]

nm : 0-0.1 [1/s]

hw : 1000-5000 [m]

cw : 0.001-0.02 [kg/m^3]

rainfall_frequency : 0.1 - 24 [number of storms of 1 hour duration per day]

Private functions#

GridProcess._regInterp(field)[source]#

Perform bilinear interpolation of field on the regular grid to unstructured 2D mesh.

Parameters:: field – data to interpolate of size m x n
Returns:: ufield field interpolated to unstructured nodes

GridProcess._buildRegGrid()[source]#

Builds the regular grid based on nodes coordinates and instantiates two interpolation objects.

The first one uses SciPy cKDTree to interpolate values from the unstructured mesh onto the regular grid based on an inverse weighting distance approach.
The second performs a bilinear interpolation from the regular grid to unstructured 2D mesh.

Note

Here the KDTree is not kept in memory, instead we store the interpolation information, namely the indices of the neighbouring nodes, and the weights of each node in the neighborhood (based on the distance). This implies that the initial distribution of the mesh coordinates remains fixed over the simulation time window.

GridProcess._updateTe()[source]#: Finds current elastic thickness map for the considered time interval.

GridProcess._cptFlex2D(dZ)[source]#: Compute the flexural response for 2D cases.

GridProcess._buildDHGrid()[source]#

Build the Driscoll-Healy (DH2, shape N x 2N) target grid for global flexure and precompute the interpolation weights linking the unstructured spherical mesh (self.mCoords) and the DH grid.

Forward direction (mesh -> DH): 3-D inverse-distance weighting via scipy.spatial.cKDTree on the existing Cartesian coordinates, which avoids polar / dateline singularities.

Backward direction (DH -> mesh): bilinear interpolation with precomputed integer indices and fractional offsets. The DH grid is implicitly extended at call time with a south-pole row and a wrap-around longitude column so every mesh point sits inside a 2 x 2 stencil.

Because the unstructured mesh is fixed for a given simulation, weights are built once and reused for every flexure step.

GridProcess._unstr2dh(field)[source]#: Inverse-distance interpolation of a 1-D unstructured field defined at self.mCoords onto the DH2 grid (shape (dh_N, 2*dh_N)).

GridProcess._dh2unstr(field_dh)[source]#: Bilinear interpolation of a DH2-grid array back to the unstructured mesh. The grid is padded at call time with a south-pole row (mean of the southernmost DH row) and a wrap-around longitude column.

GridProcess._cmptFlexGlobal(erodep_dh, te_dh, rho_infill=0.0, max_iter=50, flex_tol=0.0005, relax=1.0, anderson_depth=5)[source]#

Spectral thin-elastic-shell flexure solve on the precomputed Driscoll-Healy (DH2) grid.

Parameters:

erodep_dhnp.ndarray, shape (dh_N, 2*dh_N): Erosion(-) / deposition(+) thickness in metres on the DH2 grid. Produced by _unstr2dh() from the unstructured mesh field.
te_dhfloat or np.ndarray of shape (dh_N, 2*dh_N): Elastic thickness in metres. Scalar -> constant-Te spectral solve. Array (same grid as erodep_dh) -> iterative varying-Te solve.
rho_infillfloat: Density (kg/m^3) of the material filling the flexural moat (0 = air, 1030 = sea water).
max_iter, flex_tol, relax, anderson_depthint, float, float, int: Picard / Anderson iteration controls for the varying-Te branch. relax < 1 damps the update (useful for strong Te contrasts).

Returns:

np.ndarray, shape (dh_N, 2*dh_N): Flexural deflection in metres on the DH2 grid. Sign: negative = down (subsidence under deposition; rebound gives positive w under erosion).