landiff - a simple parallel multi-lithology non-linear diffusion model of hillslope dynamic

The problem

The simple creep diffusion equation for single lithology has the basic form:

where h(x, y, t) is the height above some arbitrary horizontal surface in the x-y plane, t is time and (x,y,t) is the transport or diffusion coefficient giving the effectiveness of the diffusion.

Using a parameter for the fraction of a given sediment, Rivenaes added a second equation to calculate the ratio of two sediments in a layer of deposited material. The modified equations include s(x, y, t) and 1-s(x,y,t) as the fraction of the two sediments.

In particular, we consider here a sedimentation scenario of n lithology. The following set of nonlinear PDEs (n+1 equations in total), derived from Rivenaes, constitute the mathematical model:

and for each sediment k in [ 1,...,n-1]:

In the above model, denote the diffusion coefficients for sediment k. In addition, is the compaction ratio of sediment-k type. Moreover, A is a constant representing the thickness of a prescribed top layer, in which sediments are transported.

The initial conditions are of the form and . As boundary conditions, most of the boundary has the no-flow condition, i.e., the homogeneous Neumann boundary condition:

Temporal discretisation

We note the time step size. Let superscript l be the time level index, such that denotes and denotes . Then, the temporal derivatives are simply approximated as

The remaining task of temporal discretisation is to choose time level l or l+1, or a combination of both, at which the right-hand-side terms of the PDEs system are to be evaluated. Different strategies will give rise to fully-explicit, semi-implicit and fully-implicit schemes. Here we compute at the fully-explicit scheme.

Fully-explicit scheme

To avoid solving systems of nonlinear algebraic equations, the right-hand-side terms of the PDEs system can use the already computed h and values. More specifically, the equations are transformed as follows, by a fully-explicit temporal discretisation:

and for each sediment k in [1,...,n-1]:

It should be noted that h is to be updated before s during each time step. This is why the newly computed (from the second equation) is immediately used to compute . Another remark is that , instead of , is used in the term on the left-hand side of the last equation. Numerical experiments show that this simple trick improves the numerical stability of this fully-explicit scheme, in which both and are computed straightforwardly. The scheme has first-order accuracy in time.

Spatial discretisation

We choose finite differences to carry out the spatial discretisation. This is mostly motivated by the numerical and programming simplicity. It can be mentioned that other spatial discretisation techniques, such as finite elements, can also use the same temporal discretisation discussed above.

Treatment of diffusion

It is standard to use centred difference for the two diffusion terms on the right-hand side of equation 2, for obtaining second-order accuracy in space. For example, centred difference applied to the term gives the following discretised form:

where the subscripts i,j are the mesh point index for a 2D uniform grid with mesh spacing and . In the above formula, the half-indexed terms are to be evaluated as, e.g., .

Treatment of convection

The equation 3 is a convection equation with respect to , because of the term. For the sake of numerical stability, one-sided upwind finite difference is preferred over centred difference, despite its first-order spatial accuracy.

To this end, it is customary to move the convection term to the left-hand side of equation when checking the flow direction. That is, gives the convection velocity. The x-component, , is approximated by , the sign of which determines how the x-component of the convection term is discretised by one-sided upwind difference. More specifically, the

term is approximated by

if we have . Otherwise, the following approximation is used:

Treatment of boundary conditions

Second-order accurate treatment of the homogeneous Neumann condition follows the standard approach by using one layer of ghost boundary points.

Numerical scheme

The complete numerical discretisation of the equation 2 in the PDEs is as follows:

where as previously defined, the half-index subscripts in the above formula mean some form of averaging for a quantity in the middle of two spatial mesh points. We typically use an arithmetic mean as follows:

For lithology k, equation 3 the complete numerical discretisation reads:

with condition (c1) corresponding to and (c2) otherwise. Similarly condition (c3) is used in case and (c4) otherwise.

Parallelisation

The parallelisation is simply done by splitting the regular grid between processors using either the rows or columns.

Examples

Two examples are provided for testing purposes. You will need to manually edit the main files to change some specific variables such as the number of lithology or the coefficient of diffusion for marine and aerial environment.

The results are CSV files that could be open in Paraview... or anything else!

References

• S.R. Clark, W. Wei, X. Cai, 2010. Numerical analysis of a dual-sediment transport model applied to Lake Okeechobee, Florida. in: Proceedings of the 9th International Symposium on Parallel and Distributed Computing, IEEE Computer Society Press, pp. 189-194.

• J.C. Rivenaes, 1992. A computer simulation model for siliclastic basin stratigraphy. Ph.D. thesis, University of Trondheim.

• J.C. Rivenaes, 1997. Application of a dual-lithology, depth-dependent diffusion equation in stratigraphic simulation. Basin Research, 4 (2), 133-146.

• W. Wei, S.R. Clark, H. Su, M. Wen, X. Cai, 2012. Balancing efficiency and accuracy for sediment transport simulations, link.

• H. Su, N. Wu, M. Wen, C, Zhang, X. Cai, 2013. Performance of Sediment Transport Simulations on NVIDIA's Kepler Architecture, in: International Conference on Computational Science, ICCS, Procedia Computer Science 18, 1275--1281.