Class FAMesh#

class flow.flowplex.FAMesh(*args, **kwargs)[source]#

This class calculates drainage area in an implicit, iterative manner using PETSc solvers. It accounts for multiple flow direction paths (SFD to MFD) based on user input declaration.

Note

The class follows the parallel approach described in Richardson et al., 2014 where the iterative nature of the computational algorithms used to solve the linear system creates the possibility of accelerating the solution by providing an initial guess.

For drainage computation, the class uses a depression-less surface and computes river incision expressed using a stream power formulation function of river discharge and slope.

If the user has turned-on the sedimentation capability, this class solves implicitly the stream power formulation accounting for a sediment transport/deposition term (Yuan et al, 2019).

Methods

`erodepSPL`()	Modified stream power law model used to represent erosion by rivers also taking into account the role played by sediment in modulating erosion and deposition rate.
`flowAccumulation`()	This function is the main entry point for flow accumulation computation.
`matrixFlow`(flowdir[, dep])	This function defines the flow direction matrices.

Initialise

__init__(*args, **kwargs)

The initialisation of FAMesh class consists in the declaration of PETSc vectors and matrices.

Public Methods

`erodepSPL`()	Modified stream power law model used to represent erosion by rivers also taking into account the role played by sediment in modulating erosion and deposition rate.
`flowAccumulation`()	This function is the main entry point for flow accumulation computation.
`matrixFlow`(flowdir[, dep])	This function defines the flow direction matrices.

Private Methods

`_buildFlowDirection`(h[, down])	This function builds from neighbouring slopes the flow directions.
`_coupledEDSystem`(eMat)	Setup matrix for the coupled linear system in which the SPL model takes into account sediment deposition.
`_distributeDownstream`(pitVol, FA, hl, step)	In cases where rivers flow in depressions, they might fill the sink completely and overspill or remain within the depression, forming a lake.
`_eroMats`(hOldArray)	Builds the erosion matrices used to solve implicitly the stream power equations for the river and ice processes.
`_getEroDepRate`()	This function computes erosion deposition rates in metres per year.
`_matrix_build`([nnz])	Creates a sparse PETSc matrix.
`_matrix_build_diag`(V[, nnz])	Builds a PETSc diagonal matrix based on a given array V
`_solve_KSP`(guess, matrix, vector1, vector2)	PETSc scalable linear equations solvers (KSP) component provides Krylov subspace iterative method and a preconditioner.

Public functions#

FAMesh.erodepSPL()[source]#

Modified stream power law model used to represent erosion by rivers also taking into account the role played by sediment in modulating erosion and deposition rate.

It calls the private function _getEroDepRate described above. Once erosion/deposition rates have been calculated, the function computes local thicknesses for the considered time step and update local elevation and cumulative erosion, deposition values.

FAMesh.flowAccumulation()[source]#

This function is the main entry point for flow accumulation computation.

Note

Flow accumulation (FA) calculations are a core component of landscape evolution models as they are often used as proxy to estimate flow discharge, sediment load, river width, bedrock erosion, and sediment deposition. Until recently, conventional FA algorithms were serial and limited to small spatial problems.

goSPL model computes the flow discharge from FA and the net precipitation rate using a parallel implicit drainage area (IDA) method proposed by Richardson et al., 2014 but adapted to unstructured grids.

It calls the following private functions:

_buildFlowDirection
_solve_KSP
_distributeDownstream

FAMesh.matrixFlow(flowdir, dep=None)[source]#

This function defines the flow direction matrices.

Note

The matrix is built incrementally looping through the number of flow direction paths defined by the user. It proceeds by assembling a local Compressed Sparse Row (CSR) matrix to a global PETSc matrix.

When setting up the flow matrix in PETSc, we preallocate the non-zero entries of the matrix before starting filling in the values. Using PETSc sparse matrix storage scheme has the advantage that matrix-vector multiplication is extremely fast.

The matrix coefficients consist of weights (comprised between 0 and 1) calculated based on the number of downslope neighbours and proportional to the slope.

Parameters:

flow – boolean to compute matrix for either downstream water or sediment transport
dep – deposition flux coefficient in case where the sediment transport/deposition term is considered.

Private functions#

FAMesh._buildFlowDirection(h, down=True)[source]#

This function builds from neighbouring slopes the flow directions. It calls a fortran subroutine that locally computes for each vertice:

the indices of receivers (downstream) nodes depending on the desired number of flow directions (SFD to MFD).
the distances to the receivers based on mesh resolution.
the associated weights calculated based on the number of receivers and proportional to the slope.

Parameters:

h – elevation numpy array
down – boolean to indicate whether the filled elevation needs to be considered or not.

FAMesh._coupledEDSystem(eMat)[source]#

Setup matrix for the coupled linear system in which the SPL model takes into account sediment deposition.

Note

The approach follows Yuan et al, 2019, where the deposition flux depends on a deposition coefficient \(G\) and is proportional to the ratio between cell area \(A\) and flow accumulation \(FA\).

The approach considers the local balance between erosion and deposition and is based on sediment flux resulting from net upstream erosion.

\[\mathrm{\frac{\eta_i^{t+\Delta t}-\eta_i^t}{\Delta t}} = \mathrm{-\kappa P^d_i \sqrt{Q_i} \frac{\eta_i^{t+\Delta t} - \eta_{rcv}^{t+\Delta t}}{\lambda_{i,rcv}}} + \mathrm{G' Q_{s_i} / \Omega_i}\]

where \(\mathrm{\lambda_{i,rcv}}\) is the length of the edges connecting the considered vertex to its receiver and \(\mathrm{\Omega_i}\) is the area (voronoi) of the node \(i\).

\(\mathrm{Q_{s_i}}\) is the upstream incoming sediment flux in m3/yr and \(\mathrm{G'}\) is equal to \(\mathrm{G \Omega_i / \bar{P}A}\).

The upstream incoming sediment flux is obtained from the total sediment flux \(\mathrm{Q_{t_i}}\) where:

\[\mathrm{Q_{t_i}^{t+\Delta t} - \sum_{ups} w_{i,j} Q_{t_u}^{t+\Delta t}}= \mathrm{(\eta_i^{t} - \eta_i^{t+\Delta t}) \frac{\Delta t}{\Omega_i}}\]

which gives:

\[\mathrm{Q_{s_i}} = \mathrm{Q_{t_i}} - \mathrm{(\eta_i^{t} - \eta_i^{t+\Delta t}) \frac{\Delta t}{\Omega_i}}\]

This system of coupled equations is solved implicitly using PETSc by assembling the matrix and vectors using the nested submatrix and subvectors and by using the fieldsplit preconditioner combining two separate preconditioners for the collections of variables.

Parameters:: eMat – erosion matrix (from the simple SPL model)

FAMesh._distributeDownstream(pitVol, FA, hl, step, ice=False)[source]#

In cases where rivers flow in depressions, they might fill the sink completely and overspill or remain within the depression, forming a lake. This function computes the excess of water (if any) able to flow dowstream.

Important

The excess water is then added to the downstream flow accumulation (FA) and used to estimate rivers’ erosion.

Parameters:

pitVol – volume of depressions
FA – excess flow accumulation array
hl – current elevation array
step – downstream distribution step
ice – boolean indicating where the ice flow is considered or not.

Returns:

pitVol, excess, nFA (updated volume in each depression, boolean set to True is excess flow remains to be distributed and new flow accumulation values)

FAMesh._eroMats(hOldArray)[source]#

Builds the erosion matrices used to solve implicitly the stream power equations for the river and ice processes.

Parameters:: hOldArray – local elevation array from previous time step
Returns:: eMat, PA where the first is a sparse PETSc matrices related to river and glacial erosion and PA is the accumulation rate.

FAMesh._getEroDepRate()[source]#

This function computes erosion deposition rates in metres per year. This is done on the filled elevation. We use the filled-limited elevation to ensure that erosion/deposition is not going to be underestimated by small depressions which are likely to be filled (either by sediments or water) during a single time step.

The simplest law to simulate fluvial incision is based on the detachment-limited stream power law, in which erosion rate depends on drainage area \(A\), net precipitation \(P\) and local slope \(S\) and takes the form:

\[E = − \kappa P^d (PA)^m S^n\]

\(\kappa\) is a dimensional coefficient describing the erodibility of the channel bed as a function of rock strength, bed roughness and climate, \(d\), \(m\) and \(n\) are dimensionless positive constants.

A similar approach is used to compute ice induced erosion where the ice flow accumulation is defined based on downstream nodes and is smoothed to better represent the erosion induced by glaciers. The ice-induced erosion uses the stream power law equation with a eordibility coefficient which is user defined. Under glacier terminus point, melted glacier flow is added to river flow accumulation likewise is the glacier-induced transported sediment flux.

Default formulation assumes \(d = 0\), \(m = 0.5\) and \(n = 1\). The precipitation exponent \(d\) allows for representation of climate-dependent chemical weathering of river bed across non-uniform rainfall.

Important

In goSPL, the coefficient n is fixed and the only variables that the user can tune are the coefficients m, d and the erodibility \(\kappa\).

The erosion rate is solved by an implicit time integration method, the matrix system is based on the receiver distributions and is assembled from local Compressed Sparse Row (CSR) matrices into a global PETSc matrix. The PETSc scalable linear equations solvers (KSP) is used with both an iterative method and a preconditioner and erosion rate solution is obtained using PETSc Richardson solver (richardson) with block Jacobian preconditioning (bjacobi).

An alternative method to the detachment-limited approach proposed above consists in accounting for the role played by sediment in modulating erosion and deposition rates. It follows the model of Yuan et al, 2019, whereby the deposition flux depends on a deposition coefficient \(G\) and is proportional to the ratio between cell area \(\mathrm{\Omega}\) and water discharge \(\mathrm{Q}=\bar{P}A\).

FAMesh._matrix_build(nnz=(1, 1))[source]#

Creates a sparse PETSc matrix.

Note

To achieve good performance during matrix assembly, the function preallocates the matrix storage by setting the array nnz.

Parameters:: nnz – array containing the number of nonzeros in the various rows
Returns:: sparse PETSc matrix

FAMesh._matrix_build_diag(V, nnz=(1, 1))[source]#

Builds a PETSc diagonal matrix based on a given array V

Parameters:

V – diagonal data array
nnz – array containing the number of nonzero blocks

Returns:

sparse PETSc matrix

FAMesh._solve_KSP(guess, matrix, vector1, vector2)[source]#

PETSc scalable linear equations solvers (KSP) component provides Krylov subspace iterative method and a preconditioner. Here, flow accumulation solution is obtained using PETSc Richardson solver (richardson) with block Jacobian preconditioning (bjacobi).

Note

The solver choice was made based on the convergence results from Richardson et al. (2014) but can be changed if better solver and preconditioner combinations are found.

Using such iterative method allows for an initial guess to be provided. When this initial guess is close to the solution, the number of iterations required for convergence dramatically decreases. Here the flow discharge solution from previous time step can be passed as an initial guess to the solver as discharge often exhibits little change between successive time intervals.

Parameters:

guess – Boolean specifying if the iterative KSP solver initial guess is nonzero (when provided it corresponds to the previous flow discharge values)
matrix – PETSc sparse matrix used by the KSP solver composed of diagonal terms set to unity (identity matrix) and off-diagonal terms (weights between 0 and 1). The weights are calculated based on the number of downslope neighbours (based on the chosen number of flow direction directions) and are proportional to the slope.
vector1 – PETSc vector corresponding to the local volume of water available for runoff during a given time step (e.g. voronoi area times local precipitation rate)
vector2 – PETSc vector corresponding to the unknown flow discharge values

Returns:

vector2 PETSc vector of the new flow discharge values