Eroders#

“Eroder” classes implement various erosion processes. They all provide an erode method for computing erosion and/or updating the elevation of the topographic surface during one time step.

Channel erosion#

Defined in fastscapelib/eroders/spl.hpp.

template<class FG, class S = typename FG::xt_selector>
class spl_eroder#

Bedrock channel erosion modelled using the Stream Power Law.

It numerically solves the Stream Power Law (SPL) using an implicit finite difference scheme 1st order in space and time. The method is detailed in Braun and Willet’s (2013) and has been slightly adapted.

SPL is defined as:

\[ \frac{\partial h}{\partial t} = - K A^m (\nabla h)^n \]

where \(K\) is an erosion coefficient, \(A\) is the upslope contributing (or drainage) area and \(\nabla h\) is the local gradient of the topographic surface.

For the linear case \(n = 1\), solving the equation is trivial. For the non-linear case \(n \neq 1\), the Newton-Raphson method is used to find the optimal solution.

Solving the SPL equation requires setting boundary conditions. In this implementation, the erosion at the base level nodes of the flow graph is set to zero.

See Braun and Willett [BW13] for more details on the numerical scheme used here.

Note

This implementation supports multiple direction flow, although only for the linear case.

Note

This implementation also prevents the formation of new closed depressions (lakes). The erosion within existing closed depressions of the input topographic surface will depend on the given flow graph:

  • If the flow graph handles routing accross those depressions, no erosion will take place within them (lakes).

  • If the flow graph does not resolve those depressions (no sink resolver applied), they will be eroded like the rest of the topographic surface (no lake).

Warning

The numerical scheme used here is stable but not implicit with respect to upslope contributing area (A). Using large time steps may still have a great impact on the solution, especially on transient states.

Template Parameters:

  • FG – The flow graph type.

  • S – The xtensor container selector for data array members.

Public Types

using flow_graph_type = FG#
using xt_selector = S#
using size_type = typename flow_graph_type::size_type#
using shape_type = typename flow_graph_type::shape_type#
using data_type = typename flow_graph_type::data_type#
using data_array_type = xt_array_t<xt_selector, data_type>#

Public Functions

template<class K>
inline spl_eroder(FG &flow_graph, K &&k_coef, double area_exp, double slope_exp, double tolerance = 1e-3)#

Create a new SPL eroder.

Parameters:

  • flow_graph – The flow graph instance

  • k_coef – The \(K\) value (spatially uniform or variable)

  • area_exp – The \(A\) exponent value

  • slope_exp – The \(\nabla h\) exponent value

  • tolerance – The tolerance of Newton-Raphson convergence for the non-linear case.

Template Parameters:

K – Either a scalar float type or an array (xtensor expression) type.

inline const data_array_type &k_coef()#

Return the \(K\) value.

template<class T>
inline void set_k_coef(T value, typename std::enable_if_t<std::is_floating_point<T>::value>* = 0)#

Set a spatially uniform, scalar \(K\) value.

template<class K>
inline void set_k_coef(K &&value, typename std::enable_if_t<xt::is_xexpression<K>::value>* = 0)#

Set a spatially variable, array \(K\) value.

The array expression or container must have the same shape than the grid arrays.

inline double area_exp()#

Return the \(A\) exponent value.

inline void set_area_exp(double value)#

Set the \(A\) exponent value.

inline double slope_exp()#

Return the \(\nabla h\) exponent value.

inline void set_slope_exp(double value)#

Set the \(\nabla h\) exponent value.

inline double tolerance()#

Return the tolerance controlling the convergence of the Newton-Raphson iterations for the non-linear case.

inline size_type n_corr()#

Returns the number of nodes for which erosion has been arbitrarily limited during the last computed time-step.

To ensure numerical stability, channel erosion may not lower the elevation of a node such that it reverts the slope with any of its direct neighbors. This prevents the formation of new closed depressions.

const data_array_type &erode(const data_array_type &elevation, const data_array_type &drainage_area, double dt)#

Solve SPL for one time step.

Parameters:

  • elevation – The elevation of surface topography at each grid node.

  • drainage_area – The upslope drainage area at each grid node.

  • dt – The duration of the time step.

template<class FG, class K>
spl_eroder<FG> fastscapelib::make_spl_eroder(FG &flow_graph, K &&k_coef, double area_exp, double slope_exp, double tolerance)#

Helper to create a new SPL eroder.

Parameters:

  • flow_graph – The flow graph instance

  • k_coef – The \(K\) value (spatially uniform or variable)

  • area_exp – The \(A\) exponent value

  • slope_exp – The \(\nabla h\) exponent value

  • tolerance – The tolerance of Newton-Raphson convergence for the non-linear case.

Template Parameters:

  • FG – The flow graph type.

  • K – Either a scalar float type or an array (xtensor expression) type.

Hillslope erosion#

Defined in fastscapelib/eroders/diffusion_adi.hpp.

template<class G, class S = typename G::xt_selector>
class diffusion_adi_eroder#

Hillslope erosion using linear diffusion.

It numerically solves the diffusion equation using an Alternating Direction Implicit (ADI) scheme.

The equation is given by:

\[ \frac{\partial h}{\partial t} = K \nabla^2 h \]

where \(K\) is an erosion coefficient (diffusivity) and \(\nabla^2 h\) is the local curvature of the topographic surface.

This equation implies that the amount of sediment eroded is linearly proportional to the local gradient of the topographic surface.

Note

This eroder assumes Dirichlet boundary conditions at the border nodes of the raster grid.

Warning

Only raster grids are supported.

Template Parameters:

  • FG – The raster grid type.

  • S – The xtensor container selector for data array members.

Public Types

using grid_type = G#
using xt_selector = S#
using data_type = typename grid_type::grid_data_type#
using data_array_type = xt_array_t<xt_selector, data_type>#
using data_tensor_type = xt_tensor_t<xt_selector, data_type, grid_type::xt_ndims()>#

Public Functions

template<class K>
inline diffusion_adi_eroder(G &grid, K &&k_coef, typename std::enable_if_t<is_raster_grid<G>::value>* = 0)#

Create a new diffusion ADI eroder.

Parameters:

  • grid – The input raster grid.

  • k_coef – The \(K\) value (spatially uniform or variable).

Template Parameters:

K – Either a scalar float type or an array (xtensor expression) type.

inline data_array_type k_coef()#

Return the \(K\) value.

template<class T>
inline void set_k_coef(T value, typename std::enable_if_t<std::is_floating_point<T>::value>* = 0)#

Set a spatially uniform, scalar \(K\) value.

template<class K>
inline void set_k_coef(K &&value, typename std::enable_if_t<xt::is_xexpression<K>::value>* = 0)#

Set a spatially variable, array \(K\) value.

The array expression or container must have the same shape than the raster grid arrays.

const data_array_type &erode(const data_array_type &elevation, double dt)#

Solve diffusion for one time step.

Parameters:

  • elevation – The elevation of surface topography at each grid node.

  • dt – The duration of the time step.

template<class G, class K>
diffusion_adi_eroder<G> fastscapelib::make_diffusion_adi_eroder(G &grid, K &&k_coef)#

Helper to create a new diffusion ADI eroder.

Parameters:

  • grid – The input raster grid.

  • k_coef – The \(K\) value (spatially uniform or variable).

Template Parameters:

  • G – The raster grid type.

  • K – Either a scalar float type or an array (xtensor expression) type.