advertisement

WATER RESOURCES RESEARCH, VOL. 38, NO. 7, 10.1029/2001WR000721, 2002 An implicit finite difference method for drainage basin evolution Sergio Fagherazzi,1 Alan D. Howard, and Patricia L. Wiberg Department of Environmental Sciences, University of Virginia, Charlottesville, Virginia, USA Received 15 June 2001; revised 25 October 2001; accepted 25 October 2001; published 24 July 2002. [1] In drainage basin evolution models the implementation of sediment transport and sediment balance is constrained by numerical instabilities. In order to overcome these difficulties an implicit finite difference method is proposed. The method simultaneously adjusts the elevation at each node of the numerical grid, is unconditionally stable, and significantly reduces the computational time. The performance of the method is tested and INDEX TERMS: 1824 Hydrology: Geomorphology briefly discussed with a numerical example. (1625); 1815 Hydrology: Erosion and sedimentation; 3230 Mathematical Geophysics: Numerical solutions; 3210 Mathematical Geophysics: Modeling; KEYWORDS: drainage basin, landscape evolution, sediment transport, detached limited model, river network 1. Introduction [2] In recent years, computational models based on partial differential equations have been developed to study the evolution of drainage basin landscapes. These models mimic key processes of landscape evolution, such as sediment transport in alluvial and bedrock channels, soil creep, rain splash, with a few simplified equations derived from general conservation principles. The equations are discretized on numerical grids and suitable algorithms link the quantities at each node. [3] In the simplest models, only two processes are present: a diffusive soil creep plus advective sediment transport in channels [Willgoose et al., 1991a; Rinaldo et al., 1995; Howard, 1994]. The alluvial sediment transport is formulated as a function of drainage area, a surrogate for discharge, and slope. When calculating sediment budgets at each node, the high variability of drainage area amplifies small variations in bottom slope leading to a solution that requires very small time steps for numerical stability, thus limiting the model applicability for large spatial and temporal domains. [4] Different techniques have been developed to overcome this obstacle, but none of them have the flexibility and the rigor that such an important problem requires. Rinaldo et al. [1995] modeled the alluvial part with a threshold-limiting formulation, in which sediment is removed in the landscape until the bed stress is less than the critical shear stress, but deposition does not occur in the alluvial channels, and the model cannot follow the fate of sediment once it arrives in the channel network. Willgoose et al. [1991] modified the equations utilizing a local asymptotic solution, but the improvement is limited, and the modified equation depends only on the elevation in the neighborhood of the considered node, whereas all the nodes of an alluvial channel have to be modified synchronously to have a realistic solution. [5] In this direction the approach of Howard [1994] is more sound. A steady state solution is utilized to adjust at the same time the slope of all the nodes belonging to a channel. However, the method involves extensive bookkeeping, and since it refers to steady state channel gradients, it is unable to follow transient dynamics. [6] Here a new implicit finite difference method is presented. The method solves the original equations for the alluvial transport adjusting simultaneously the bottom elevation in each node of the domain. 2. Method [7] In the present analysis we focus on the sediment transport and drainage basin evolution due to alluvial channels. This notwithstanding, the same formulation can be applied, with minor modifications, to other processes responsible for landscape form and evolution. [8] In drainage basin models the basin area is partitioned in elements, with a node at the center of each element. Variations in elevation are then computed with a mass balance in each node: @z* 1 X ¼ Q*sb ; @t* a* where z* is the elevation, a* is the basin area pertaining to the node, Q*sb is the volumetric sediment discharge entering or leaving the node from one of the neighboring nodes (the superscript asterisk indicates dimensional quantities). The sediment discharge can be derived from sediment transport laws and expressed in the form [Howard, 1994; Willgoose et al., 1991a] Q*sb ¼ KA*m S*p ¼ KA*m 1 Now at Department of Geological Sciences and School of Computational Science and Informational Technology, Florida State University, Tallahassee, Florida, USA. Copyright 2002 by the American Geophysical Union. 0043-1397/02/2001WR000721$09.00 ð1Þ z* p ; x* ð2Þ where K is a sediment transport parameter, A* is the drainage area at the point, S* is the bottom slope, p and m are constant coefficients, x is the distance between two nodes. The area A* is a surrogate variable for the mean annual discharge at a site and it is determined in the 21 - 1 21 - 2 FAGHERAZZI ET AL.: TECHNICAL NOTE the method can be easily extended to models where discharge and sediments is coming out in more than one direction [Tarboton, 1997; Costacabral and Burges, 1994]. [10] In the standard, explicit scheme used to solve equation (4) the right term is evaluated at the current time n, while the time derivative is expressed as a forward difference in time. Naming the right-hand term of equation (4) f, the explicit scheme is znþ1 zni i ¼ f ðzn ; An Þ; t where the quantities calculated at the time step n are utilized to calculate the elevation at the time step n + 1. Howard [1994] and Willgoose et al. [1991a, 1991b] pointed out that this method requires very small time steps to be stable, with a high computational burden. Here we explore the performance of an implicit method, evaluating the righthand term at time n + 1. A first simplification adopted consists in calculating the drainage area at time n, in that way only the slope in equation (4) is an unknown function of the elevation at time n + 1. Moreover, we linearize the function f with a Taylor expansion around its value at t = n. The implicit method can then be written as Figure 1. Computational mesh. following way: first, each node is linked to one of the neighboring nodes following the steepest gradient; the node is thus draining in this direction. Second the drainage area is computed as the sum of the areas pertaining to nodes that, directly or indirectly, drain in the considered one. The alluvial sediment is then supposed to be transported following the drainage directions. For simplicity, in equation (2) we did not introduce a threshold shear stress for sediment transport, but the method presented is valid for any formulation of the sediment transport law. As by Willgoose et al. [1991b], we nondimensionalize equations (1) and (2) utilizing the following expressions: x* ¼ Lx; a* ¼ L2 a; Z A* ¼ L2 A; S* ¼ S; t* ¼ Tt L z* ¼ Zz; ð3Þ where z, x, a, A, S, and t are the nondimensional elevation, spatial coordinate, basin area pertaining to the node, drainage area, slope, and time, respectively. Z, L, and T are the corresponding scales for elevation, length, and time. [9] Equation (1) can then be expressed as a function of elevation and drainage area of the neighboring points: @zi X H Am zj zi p H Am zi zk p j i ¼ @t a xj a xi j ð4Þ with the nondimensional number H equal to H ¼ KZ p1 TL2mp2 : ð6Þ ð5Þ In equation (4) j are the points that are draining to the point i, and k is the point where i is draining to (Figure 1). The number j depends on how the drainage area is aggregated on the grid. In the formulation adopted here we have an outgoing sediment discharge to only one of the neighboring nodes, but X @f znþ1 zni i nþ1 n n n znþ1 zn ; ð7Þ ¼ f ðz ; A Þ f ðz ; A Þ þ j j @zj n t j applying this formulation to equation (4) and rearranging the terms, we derive an equation for each node i, " 1 þ t X !# Fjn þ Fin t znþ1 i j X Fjn znþ1 tFin znþ1 ¼ zni j k j þtC þ t " X # Fjn zni znj þ Fin zni znk ; ð8Þ j where the quantities Fjn, Fin, and C are evaluated at the time step n and are Fjn ¼ n zj zni p1 @f pHAm j ¼ xj @zj n axj Fin ¼ C¼ n @f pHAm zi znk p1 i ¼ @zk n axi xi ð9Þ n n X HAm zj zni p HAm zi znk p j i : a xj a xi j If R is the number of nodes, equation (8) represents a system of R linear equations in the unknown elevations zin + 1 at time n + 1. The R R matrix related to the system is symmetric and definite positive. Furthermore, the matrix is sparse since each node is linked only to the neighboring nodes, so not more than nine elements among the R elements in the line i of the matrix are different than zero. The storing of the matrix is effectively performed by retaining only the nonzero elements and their location. The system is then solved with a preconditioned conjugate FAGHERAZZI ET AL.: TECHNICAL NOTE 21 - 3 Figure 2. Initial and boundary conditions for the test case. gradient method, where the preconditioning matrix is calculated with the Kershaw’s incomplete factorization of the original matrix decomposed with the Cholesky method [Axelsson, 1994]. 3. Results [11] The method is unconditionally stable; however, errors are introduced when a large time step is adopted. As a matter of fact, the linearization of the right-hand term in equation (7) is only valid for small time steps, when changes in elevation are limited. Most importantly, the method utilizes a ‘‘frozen’’ drainage area each time step, without tracking its variations caused by elevation changes. To test the validity of the method, we compare results obtained with the implicit and the explicit schemes for a drainage basin evolution example. The domain is a 100 100 square Cartesian grid, and the length scale L is equal to the mesh size. The parameters utilized in the simulation are H = 100, m = 1.4, and p = 2.1, corresponding to the parameters by Howard [1994, Figure 18]. The starting topography is a tilted plane in the vertical direction with a constant gradient of 0.05 to which is added small random elevations selected from a normal distribution with zero mean and standard deviation equal to 0.1. The upper, left, and right boundaries are impermeable (neither sediment flux Figure 3. Number of nodes changing drainage direction in each time step. The data are averaged over 20 steps to eliminate short-term fluctuations. Figure 4. Simulation 1: (a) initial conditions, (b) changes in elevation distribution at t = 0.005 using the explicit formulation (10000 iterations), and (c) changes in elevation distribution at t = 0.005 using the implicit formulation (10 iterations). 21 - 4 FAGHERAZZI ET AL.: TECHNICAL NOTE Figure 5. (a) Error of the implicit method as a function of the time step for the test case shown in Figure 2. Open circles refer to the first simulation from t = 0 to t = 0.0005 (Figure 4), while the solid circles refer to the second simulation from t = 0.0005 to t = 0.001 (Figure 6). The initial condition for the second simulation is the result of the explicit method for the first simulation. (b) ratio of the CPU time used by the implicit method relative to the CPU time used by the explicit method. nor discharge are allowed through them), whereas the lower boundary is kept at a constant elevation (Figure 2). [12] The explicit method requires a time step less than or equal to 5 107 to be stable. We run the explicit method for a duration equal to 0.01. Plotting the number of sites that change drainage direction in each time step as a function of time, we clearly see that most of the changes occur between t = 0 and t = 0.0025 (Figure 3). In this period the network develops from the random distribution of elevations. Small differences can lead to very different final results, as already pointed out by Willgoose et al. [1991a]. [13] The solution is then compared with the solution of the implicit method with the same initial conditions. In the implicit method we utilize different time steps, ranging from 0.005 (1 iteration) to 2.5 106 (2000 iterations). Assuming that the explicit method returns the exact solution, the error is calculated as the standard deviation of the implicit result with respect to the explicit one. The error is then normalized dividing it by the standard deviation of the elevations at the end of the simulation with respect to the initial elevations. Figure 6. (opposite) Simulation 2: (a) initial conditions, corresponding to the result of the first simulation with the explicit method, (b) changes in elevation distribution at t = 0.01 using the explicit formulation (10000 iterations), and (c) changes in elevation distribution at t = 0.01 using the implicit formulation (10 iterations). FAGHERAZZI ET AL.: TECHNICAL NOTE [14] The first simulation is carried out from t = 0 to t = 0.005 (Figure 3); the initial condition is a flat surface where erosion has not formed drainage patterns yet (Figure 4a). The error decreases with the reduction of the time step in the implicit method, corresponding to an increase in the number of iterations (Figure 5a, open circles). The high error value for the implicit scheme with few iterations is due to the development of a different channel network on the surface (Figures 4b and 4c). Differences arise from the fact that the discharge is kept ‘‘frozen’’ in the implicit method and does not vary during the simulation. Still the result is qualitatively similar, and the CPU time used is considerably less than the CPU time used for the explicit. For a more accurate result the number of iterations in the implicit scheme has to be increased, with consequent increment of CPU time (Figure 5b). A second simulation was performed running the two models from t = 0.005 to t = 0.01 utilizing the elevations calculated with the explicit method at time t = 0.005 as initial conditions (Figure 6a). In this way, alluvial erosion acts in a channel network already developed, and drainage area variations are almost nonexistent (Figures 6b and 6c). Now the error for the implicit method is <6% even with a single iteration (Figure 5a, solid circles). [15] In situations where the drainage area does not vary significantly, as between t = 0.005 and t = 0.01 (Figure 3), the superiority of the implicit method is clear, with an utilization of CPU time <1% of that required for the explicit formulation. The rate of change in drainage area (Figure 3) might then be utilized to modify in an adaptive way the time step of the implicit method. For example, in the first simulation we can start with a time step equal to 5 106 for 200 iterations, then switch to a time step of 5 104 for the remaining 9 iterations. The resulting error is <7%, and the computational time is one tenth of the explicit one. Other situations where the implicit method has to be carefully applied are in locations where alluvial deposition exceeds alluvial erosion. In these cases (as, for example, in alluvial fans), deposition increases local elevations, resulting in a high probability of channel avulsion. 21 - 5 models presents several advantages. The method is unconditionally stable, conserving sediment mass for any time step adopted. The formulation can be easily extended to incorporate other important processes in landscape evolution, such as soil creep and rain splash diffusion. The method can be applied both to Cartesian grids and unstructured triangular grids. Computationally, the implicit scheme is very efficient, reducing the CPU time by up to a factor 100 with respect to an explicit formulation. Since the scheme does not update the drainage area during a time step, particular attention has to be paid in situations where the distribution of drainage area changes fast. In these cases, such as the first stages of channalization of flat surfaces or channel avulsion in alluvial depositional environments, a small time step in the implicit scheme is necessary to follow the transient dynamics. Limitations in the method applicability may finally arise when modeling mixed grain sizes routing and sediment stratigraphy. [17] Acknowledgments. Support for this study was provided by the Office of Naval Research, Geoclutter Program, grant N00014-00-1-0822. References Axelsson, O., Iterative Solution Methods, 672 pp., Cambridge Univ. Press, New York, 1994. Costacabral, M. C., and S. J. Burges, Digital elevation model networks (DEMON): A model of flow over hillslopes for computation of contributing and dispersal areas, Water Resour. Res., 30(6), 1681 – 1692, 1994. Howard, A. D., A detachment-limited model of drainage-basin evolution, Water Resour. Res., 30(7), 2261 – 2285, 1994. Rinaldo, A., W. E. Dietrich, R. Rigon, G. K. Vogel, and I. Rodriguez-Iturbe, Geomorphological signatures of varying climate, Nature, 374, 632 – 635, 1995. Tarboton, D. G., A new method for the determination of flow directions and upslope areas in grid digital elevation models, Water Resour. Res., 33(2), 309 – 319, 1997. Willgoose, G., R. L. Bras, and I. Rodriguez-Iturbe, A coupled channel network growth and hillslope evolution model, 1, Theory, Water Resour. Res., 27(7), 1671 – 1684, 1991a. Willgoose, G., R. L. Bras, and I. Rodriguez-Iturbe, A coupled channel network growth and hillslope evolution model, 2, Nondimensionalization and applications, Water Resour. Res., 27(7), 1685 – 1696, 1991b. 4. Conclusions [16] The use of an implicit finite difference method for alluvial sediment transport in drainage basin evolution S. Fagherazzi, Department of Geological Sciences, Florida State University, Tallahassee, FL 32306, USA. (fagherazzi@gly.fsu.edu) A. D. Howard and P. L. Wiberg, Department of Environmental Sciences, University of Virginia, Charlottesville, VA 22903, USA.