Hillslope and marine deposition

Hillslope: soil creep

Hillslope processes are known to strongly influence catchment morphology and drainage density and several formulations of hillslope transport laws have been proposed. Most of these formulations are based on a mass conservation equation and assume that a layer of soil available for transport is always present (i.e. precluding case of bare exposed bedrock) and that dissolution and mass transport in solution can be neglected.

Under such assumptions and via the Exner’s law, the mass conservation equation widely applied in landscape modelling is of the form:

\[\mathrm{\frac{\partial \eta}{\partial t}} = -\mathrm{\nabla \cdot {q_{ds}}}\]

where \(\mathrm{q_{ds}}\) is the volumetric soil flux of transportable sediment per unit width of the land surface. In its simplest form, \(\mathrm{q_{ds}}\) obeys the Culling model and hypothesises a proportional relationship to local hillslope gradient (i.e. \(\mathrm{q_{ds}=-D\nabla \eta}\), also referred to as the creep diffusion equation):

\[\mathrm{\frac{\partial \eta}{\partial t}} = \mathrm{D \nabla^2 \eta}\]

in which \(\mathrm{D}\) is the diffusion coefficient that encapsulates a variety of processes operating on the superficial soil layer. As an example, \(\mathrm{D}\) may vary as a function of substrate, lithology, soil depth, climate and biological activity.

In gospl, hillslope processes rely on this approximation even though field evidence suggest that the creep approximation is only rarely appropriate.

For a discrete element, considering a node \(\mathrm{i}\) the implicit finite volume representation of the above equation is:

\[\mathrm{\frac{\partial \eta_i}{\partial t}} = \mathrm{\frac{\eta_i^{t+\Delta t}-\eta_i^t}{\Delta t} = D \sum_{j=1}^N \frac{ \chi_{i,j}(\eta_j^{t+\Delta t} - \eta_i^{t+\Delta t}) }{\Omega_i \lambda_{i,j}} }\]

\(\mathrm{N}\) is the number of neighbours surrounding node \(\mathrm{i}\), \(\mathrm{\Omega_i}\) is the voronoi area, \(\mathrm{\lambda_{i,j}}\) is the length of the edge connecting the considered nodes and \(\mathrm{\chi_{i,j}}\) is the length of voronoi face shared by nodes \(\mathrm{i}\) and \(\mathrm{j}\).

Applied to the entire domain, the equation above can be rewritten as a matrix system:

\[\mathrm{\mathbf Q \boldsymbol\eta^{t+\Delta t}} = \mathrm{\boldsymbol\eta^{t}}\]

where \(\mathrm{\mathbf Q}\) is sparse. The matrix terms only depend on the diffusion coefficient \(\mathrm{D}\), the grid parameters and voronoi variables (\(\mathrm{\chi_{i,j}}\), \(\mathrm{\lambda_{i,j}}\), \(\mathrm{\Omega_i}\)).

In gospl, these parameters remain fixed during a model run and therefore \(\mathrm{\mathbf Q}\) needs to be created only once at initialisation. At each iteration, hillslope induced changes in elevation \(\mathrm{\boldsymbol \eta}\) are then obtained in a similar way as for the solution of the other systems using PETSc Richardson solver and block Jacobi preconditioning.

Marine deposition

In the marine realm, a nonlinear diffusion model is used for sediment-transport by rivers. When the dual lithology is activated, gospl accounts for distinct transport coefficients for the two different grain sizes (Rivenaes, 1992).

Sediment transport is modelled through nonlinear diffusion equation. The rate of elevation change in the marine environment is governed by:

\[\mathrm{\frac{\partial \eta}{\partial t}} = \mathrm{\nabla \cdot \left( K_M(\eta) \nabla \eta \right)} + Q_{sr}\]

where \(\mathrm{K_M}\) is the marine sediment transport coefficient (m2/yr), and \(\mathrm{Q_{sr}}\) is the sediment flux coming at the river mouth. As the model progresses over time so does the shoreline position due to both offshore sedimentation and prescribed eustatic sea-level variations.

As already mentioned, the distinct transport efficiency of different grain sizes is included in our algorithm for marine sediment transport and deposition by distinct transport coefficients, \(\mathrm{K_{M1}}\) and \(\mathrm{K_{M2}}\) for coarse and fine sediments, respectively. gospl considers that these transport coefficients are uniform in space and in time. When using equation above, \(\mathrm{Q_{sr}}\) is divided in two components \(\mathrm{Q_{sr1}}\) and \(\mathrm{Q_{sr2}}\) which are the fully uncompacted flux of coarse and fine coming from the continental domain of the model, available for transport and deposition.

During a single time step, marine sediment transport is performed until all available sediment transported by rivers to the ocean have been diffused and the accumulations on these specific nodes remains below water depth or below a prescribed slope computed based on local water depth and distance to the nearest coastline.

Like for inland deposition, the coarser sediments are deposited first, followed by the finer ones. It allows for finer sediments to be deposited further and reproduce the standard behaviour observed in stratigraphic architectures. The implicit finite volume approach already presented above is implemented and the solution for elevation changes are obtained in a similar fashion as for the solution of the other matrix systems using PETSc Richardson solver and block Jacobi preconditioning.