Modelling grain-size distributions. A comparison of two models and their numerical solution

Tectonophysics 320 (2000) 347–373 www.elsevier.com/locate/tecto

Modelling grain-size distributions. A comparison of two models and their numerical solution Taco den Bezemer *, Henk Kooi, Jurjen Kranenborg, Sierd Cloetingh Faculty of Earth Sciences, Vrije Universiteit Amsterdam, De Boelelaan 1085, 1081 HV Amsterdam, The Netherlands

Abstract Three main diffusion-based models are currently used to study grain-size distributions. In this paper, two of these approaches — perfect sorting and imperfect sorting — are compared in a parameter study. First, the numerical solution of the imperfect-sorting model is extensively discussed, and numerical tests are performed. Then, the two sedimentation models are compared for a basin under varying conditions. For some of the imposed variations, predictions of both models differ markedly due to the different approach. The position of the gravel front in the perfect sorting model depends on gravel input and proximal accommodation space. The position of grain-size boundaries in the imperfect-sorting model is strongly controlled by gravel input, the position of the basin axis and the difference in diffusivities. As a result, those two models may predict gravel progradation for different situations. Both models suggest that gravel progradation should always be coupled with sedimentation rates in order to suggest an explanation of gravel progradation observed in the geological record. Simulations with the imperfect sorting model show that this criterion may also fail, showing that a unique interpretation of gravel progradation may be impossible. © 2000 Elsevier Science B.V. All rights reserved.

1. Introduction Numerical sedimentation models linking sediment transport and basin subsidence are playing an increasingly important role in interpreting how complex depositional processes interact through time to produce the sediment architectures observed in the sedimentary record. At present, numerical sedimentation models serve two main purposes. The first is in independently testing hypotheses conceived from natural examples. The second is educational in the sense that forward models can show a development in time, whereas in natural examples, only the final result can be * Corresponding author. Tel.: +31-54-649-1407; fax: +31-20-646-2457. E-mail address: [email protected] ( T. den Bezemer)

observed. Very often, these models are immediately used to model a certain geological example but, with a few exceptions, numerical tests or basic parameter tests are not provided because these are not the primary interest of most geologists. For modellers planning to modify or incorporate a certain model, however, these tests contain vital information. To this purpose, we will present both numerical and parameter tests of two different sedimentation models: a perfect-sorting model (Angevine et al. 1990) and an imperfect-sorting model (Rivenaes, 1992). Angevine et al. (1990) present basic tests of their model, but this was not done by Rivenaes (1992) for his model. His model, however, is probably more realistic and much more difficult to interpret. We will compare numerical techniques to solve the diffusion equations arising in the model formu-

0040-1951/00/$ - see front matter © 2000 Elsevier Science B.V. All rights reserved. PII: S0 0 4 0- 1 9 51 ( 0 0 ) 0 00 4 1 -X

348

T. den Bezemer et al. / Tectonophysics 320 (2000) 347–373

lation of the imperfect-sorting model (Rivenaes, 1992), and we will investigate model behaviour for a few basic controls like subsidence, gravel supply and water supply. Although primarily aimed at modellers, we believe that this paper also contains interesting ideas on the problem of interpretation of gravel progradation.

2. Geometrical and process-based models At present, many different approaches exist to simulate basin filling. Basin-filling models may be classified into geometrical and process-based models (Rivenaes, 1992). The distinction refers to two fundamentally different approaches. Geometrical models use empirical rules to mimic the response of the basin-filling system. In the absence of subsidence, many sediment-transport systems tend toward an equilibrium state, called the graded state in rivers, but this may exist in any transport system. A system will tend to a different equilibrium when the subsidence changes. Therefore, the surface topography adjusts itself until the rate of deposition and subsidence are in balance. Geometrical models (e.g. Jervey, 1988; Lawrence et al., 1990) make assumptions regarding this constant topography characteristic. Very often, depositional slope angles and critical bypass angles must be given explicitly and are kept constant during the simulation. Process-based models make use of conservation laws and flux formulations and are more realistic in that the sedimentation rate is not implicitly prescribed as in geometrical models. For instance, geometrical models respond instantaneously to an imposed model disturbance, while process models will respond with a certain time lag controlled by diffusivity (diffusion models). Thus, in processbased models, the geometry of the depositional system is calculated as part of the solution. Although process-based models will provide more detailed predictions than geometrical models (Angevine et al., 1990; Cross and Harbaugh, 1990), geometrical models may in many instances provide very useful first-order insights (e.g. Pitman, 1978). It is also important to note that even the most sophisticated process-based models, the fluid-

flow models, use semi-empirical rules for entrainment and deposition of sediment in relation to flow velocity and are consequently still of limited generality (Angevine et al., 1990).

3. Selective deposition versus abrasion The primary effect of creating topography is to induce erosion, mass transport and deposition. In alluvial systems transporting mixed size, coarse sediments, deposition selectively removes the coarsest material from the flow. The main cause for this down-slope sorting has been a matter of debate (e.g. Ferguson et al., 1996). At first abrasion (mechanical wear of rocks or particles; Krumbein, 1941, Shumm and Stevens, 1973) was thought to be the primary mechanism for down-slope sorting. This was primarily based on observations that clasts of a different composition, but of almost identical density, become finer at different rates (Shaw and Kellerhals, 1982), indicating that abrasion is dominant. However, studies on down-slope fining rates of alluvial fans showed that these rates were several orders of magnitude higher than abrasion rates estimated from laboratory experiments, suggesting that selective deposition is the most important mechanism. Although this could be partly biased, due to scaling effects of the laboratory experiments, other observations showed a positive relationship between the rate of fining and the rate of deposition (see Fig. 1), suggesting that selective deposition is the dominant down-slope sorting mechanism (Brierley and Hickin, 1985; Angevine et al., 1990).

4. Down-slope sorting models Probably the most natural way to simulate settling and entrainment of sediment of mixed grain size is to use some kind of hydrodynamic model (e.g. Tetzlaff and Harbaugh, 1989). In these models, individual fluid elements are traced, and grains are entrained/deposited, depending on fluid velocity, grain size and settling velocity. Tracing fluid elements is computationally very demanding, and this is primarily why diffusion models have been used more often in studying basin fill.

T. den Bezemer et al. / Tectonophysics 320 (2000) 347–373

Fig. 1. Grain size and sediment thickness as a function of downstream distance, x. D1 is the grain size relative to its value, D 0 at the upstream end of the basin. X1 is downstream distance relative to total basin length L . Open circles show grain-size 0 data, the dashed line shows the data trend, and the solid line shows the sediment thickness (after Paola, 1988).

Diffusion models do not trace the path of each grain in detail, but use mass conservation in combination with a bulk transport relation, which averages over spatial scales much larger than individual grains and includes the integrated effect of transport by several processes (e.g. river processes, creep, slump, small slides) and can therefore handle larger scales of space and time. A drawback of these models is that they do not have a natural way to simulate down-slope sorting. There are three approaches to circumvent this problem. The first approach is the simplest and was employed by Hardy and Waltham (1992) and Flemings and Grotzinger (1996). In this approach, the concept of a critical angle is used where a certain grainsize population bypasses, or deposits on, a slope that is calculated by diffusion. This class of models is partly geometrical in the sense that there is no conservation of mass of individual grain sizes. The second approach was employed by Angevine et al. (1990). They assumed perfect sorting, which means that only one grain-size population is deposited ( leaving finer material to be carried on by the flow) until, at a certain gridpoint within the model, the flux of that grainsize population (e.g. gravel ) is exhausted. The flux of gravel supplied to the basin, is a fraction of total flux of sediments supplied to the model. In

349

this model, the mass of individual grain sizes is preserved. As in the basic topographic diffusion model and in the previous approach, sediment transport is governed by topographic gradient, and grain size does not influence the depositional profile. It is important to mention that in the perfectsorting model, the rate of down-slope sorting is only controlled by the rate of deposition. The third approach involves imperfect sorting of two lithologies (Rivenaes, 1992). In this approach, two coupled diffusion equations are solved, each describing the mass balance of a single lithology. Sediment transport is again proportional to the local gradient and a lithology-specific diffusion coefficient but, unlike in simple topographic diffusion, also to the fraction of the lithology in a thin transport layer (boundary layer). Fractionation (selective deposition) of sediment occurs because one grain-size class is made more easily transportable than the other by assigning a higher diffusion coefficient to it. The model makes use of a transport boundary layer (see Appendix A), which makes transport of each lithology dependent on the availability of this lithology in a thin top layer. The grain-size fraction within this toplayer adjusts itself to the grain-size fraction, which is delivered from upslope and, if erosion occurs, to the grain-size fraction that has to be entrained from below the layer. This toplayer is a surrogate for (a stacked system of ) natural transport layers that ensure that what may be eroded adapts to that which is present in the subsurface. In contrast, in the imperfect-sorting model, this sorting rate is controlled by both the rate of deposition and the difference in diffusivities of gravel and sand. The last two approaches — the perfect-sorting diffusion model and the imperfect-sorting dual lithology model — are used to study grain-size distributions under varying geological conditions discussed later in this paper. First, their numerical solution is discussed in Section 5.

5. Numerical solution of the sedimentation models 5.1. Perfect-sorting model The perfect-sorting approach makes use of a single diffusion equation and assumes that sedi-

350

T. den Bezemer et al. / Tectonophysics 320 (2000) 347–373

ment transport is a function of only one transport coefficient. The standard diffusion equation is given by: ∂h ∂t

=k

∂2h

(1)

∂x2

where h is topgraphic height, t is time, k is diffusivity and x is a spatial dimension (see Appendix, for derivation). The standard diffusion equation is very well known and can be solved easily, using existing numerical techniques (Press et al., 1992). We have chosen to solve the system using an implicit Euler solution scheme (Press et al., 1992; Smith, 1992), making a typical model run (see next section) about a factor 10 faster, compared with solving the system using an explicit Euler scheme. This is mainly because the matrix resulting from this type of diffusion equation, is a positive definite tridiagonal matrix. Such a matrix does not need any row exchanges, and its lower–upper decomposition is LDLT. As a result, only 2n calculations are needed to solve the system (Strang, 1988). 5.2. Imperfect-sorting model The solution of the imperfect-sorting model proved considerably more tedious. The imperfectsorting model uses a system of coupled diffusion equations (see Appendix A for derivation): A

∂F ∂t

+F9

∂h ∂t

=K

+(1−F)

∂ G ∂x

∂h ∂t

A B A

=K

F

∂h

∂x

∂

S ∂x

−A

5.3. Explicit Euler solution scheme for the imperfect-sorting model

∂F ∂t

(1−F )

∂h ∂x

boundary conditions. Constant elevation–constant fraction boundary conditions are used. Throughout this paper, a constant thickness of the surface layer is adopted, so the model does not differentiate between thicknesses of transport layers of the hillslope (a regolith, which is moving by creep, slumps and slides) or river bed (suspension load and bed load ). Although the toplayer thickness is likely to vary (Rivenaes, 1992) in space and time, tests with a laterally varying A did not have any significant influence on the results. This thickness does, however, influence the numerical stability of the system. A small A (<1 cm) severely increases the computing time. As is shown below, this system is more difficult to solve than a single diffusion equation. The most common technique for solving diffusion type equations numerically is by using finite differences. This method approximates the derivatives by a truncated Taylor series (see Smith, 1992). Four different numerical solution techniques are compared. The easiest to use is the explicit Euler method (forward in time, first-order accuracy), whereas the three others, implicit Euler (backward in time, first order accuracy), Crank–Nicolson (backward in time, second-order accuracy) and B3 (three step backward in time, second-order accuracy) are more complicated (e.g. Press et al., 1992; Smith, 1992).

B

(2)

(ignoring compaction and assuming constant K G and K ). A is the thickness of the surface layer, F S is the fraction in the surface layer, t is the time, h is the topographic height, K and K are diffusion G S coefficients of gravel and sand respectively, and x is the spatial dimension. Note that one difference of this coupled set of equations, compared with the single diffusion Eq. (1), is that for Eq. (2), height changes may occur for a linear topographic height profile, when F varies along the profile. This equation set is supplemented by two

The system of non-linear partial differential Eq. (2) may be solved using the explicit Euler scheme by rewriting the set of equations, separating time derivatives of F and h. Summing the gravel and sand Eq. (2) gives an equation with the temporal evolution of height, which we will refer to as the height (h)-equation, despite the fact that it still contains terms in F on its right-hand side. The companion equation describing the temporal evolution of F, is referred to as the F-equation. The height equation is expressed by:

A B

A

B

∂h ∂ ∂h ∂ F +K (1−F ) . =K S G ∂t ∂x ∂x ∂x ∂x

∂h

(3)

T. den Bezemer et al. / Tectonophysics 320 (2000) 347–373

The new topographic height field, hn+1 , is obtained i by using Fn. The new fraction value, Fn+1 , is i i obtained by using either equation of (2), using the computed height field, hn+1. The main problem is i in calculating the time step necessary for numerical stability. For a simple linear diffusion equation like Eq. (1), the constraint on the time step is (e.g. Smith, 1992): Dt≤

Dx2

(4)

2K

where K is the diffusivity. It turns out that for the height equation obtained in Eq. (3), this time step is a fair approximation for numerical stability in h (note that for a system in which F=constant, the system in Eq. (1) simplifies to the simple diffusion equation), when K is assigned the highest value of K and K . Unfortunately, the maximum G S time step is not controlled by Eq. (3), but by Eq. (5). What happens in practice, when using the time step calculated in Eq. (3), is that the updated F-values reach values below 0 or above 1, which are physically impossible solutions. In order to obtain an idea of the time step necessary to prevent this, one has to analyse the Fequation. The F-equation is obtained by eliminating time derivatives of h from Eq. (2). ∂F ∂t

1−F9

=

A −

A B A A BB ∂

∂h

K F G ∂x ∂x

∂h

∂

F9

(1−F ) K A S ∂x ∂x

(5)

This F-equation is of the form: ∂F ∂t

=b

∂F ∂x

+cF+d.

(6)

The constants b, c and d are given by: b=

1−F9 A

+

F9 A

K

K

∂h G ∂x

+

A B

∂2h S ∂x2

∂h K A S ∂x F9

d=−

F9 A

c=

K

method basically states that there is a direction along which the integration of an equation of the form of Eq. (6) transforms to the integration of an ordinary differential equation. The new F value is given by:

A

1−F9 A

A B

∂2h . S ∂x2

A B

∂2h K G ∂x2

(7)

Eq. (6) may be solved using a technique called the ‘method of characteristics’ (Smith, 1992). This

B

d d F= F + ecDt− old c c

(8)

where F is the interpolated F value, belonging old to the back-shooted position, x , given by: old x =x+bDt (9) old and the time step is controlled by: cDt≤0.5.

(10)

The factor, c, contains the curvature d2h/dx2, and the time step is thus very hard to predict. An estimate of c may be obtained, using the fact that maximum slope dh/dx in the model in subsequent paragraphs does not exceed 1 and taking K =K G S (=max K ). With dh/dx<1 and a grid spacing G Dx, it follows that d2h/dx2<1/Dx, and the expression for c simplifies to: c=

K

S ADx

(11)

where Dx is the node distance. The time step, Dt, then equals: Dt=

.

351

ADx

(12) 2K S In most conditions, Eq. (12) is a too much of a restriction on Dt [compare this with the time step given by Eq. (4)]. An explicit Euler solution scheme is probably only useful in situations where the transport coefficients, K and K , are small, g s gradients in h are small, and temporal and spatial lengths are relatively short. 5.4. Implicit solution schemes The implicit Euler and B3 solution schemes are unconditionally stable for simple diffusion problems (Press et al., 1992; Smith, 1992) but computationally more demanding than explicit solution schemes. The excellent stability properties allow

352

T. den Bezemer et al. / Tectonophysics 320 (2000) 347–373

the use of very large time steps. The main drawback of the implicit Euler method is its first-order accuracy, caused by a local truncation error. The accuracy can be improved using small time steps and a finer mesh, but this obviously costs computer time. The B3 and the Crank–Nicolson method are both second-order-accurate. The Crank–Nicolson method is only conditionally stable, however (Smith, 1992), and its time step can be only roughly twice the time step required by an explicit solution scheme. We encoded all three implicit solutions schemes and compared their numerical behaviour. we will write out one of these schemes, the B3 scheme below. The other two schemes may be easily derived in a similar manner. The standard formulation for the B3 (three-time-step level ) method for a simple diffusion equation [see Eq. (1)] is (Smith, 1992): 3 (hj+1 −hj ) 1 (hj −hj−1 ) i i − i i 2 Dt 2 Dt (hj+1 −2hj+1 +hj+1 ) i i−1 = i+1 (Dx)2

(13)

where i denotes the spatial position and j the timestep level. The right-hand side of this equation is equal to the right-hand side of an implicit Euler scheme. The left-hand side employs three time levels, the old one ( j−1), the current one ( j), and the new one ( j+1), which is to be calculated.

hand side of Eq. (2)] gives: K (Fj+1 (hj+1 −hj+1 )+Fj+1 (hj+1 −hj+1 )) G i+1/2 i+1 i i−1/2 i−1 i . (Dx)2 (15) The difference between mass accumulation, M, and flow, Q, for every node, i, and each lithology expresses the residual, R, for gravel: =M −Q we=1, 2…N G,we G,we Gi and for sand:

R

(16)

=M −Q we=1, 2…N (17) S,we S,we S,i where N is the number of gridpoints. When the residual is zero for every node and each lithology, the equation set (1) is solved. The residual can be given as a vector expression as: R

, R …R , R ]T. (18) G,1 S,1 G,N G,N ( The superscript T indicates that the vector is transposed ). The solution vector is:

R=[R

x=[hj+1 , Fj+1 …hj+1 , Fj+1 ]T. (19) 1 1 N N The set of equations can be expressed as a vectorvalued function: R(x)=0.

(20)

This is a non-linear equation that may be solved using Newton iterations (e.g. Dennis and Schnabel, 1983):

5.5. B3 solution scheme imperfect-sorting model

J Dx=−R(xk) k with

(21)

Applying the B3 method to the mass accumulation term, M [ left-hand side of Eq. (2)], of the gravel equation gives:

Dx=xk+1−xk

(22)

M=

3A (Fj+1−Fj) A (Fj−Fj−1) i i − i i 2 Dt 2 Dt +

3F9 (hj+1−hj) F9 (hj−hj−1) i i − i i . 2 Dt 2 Dt

and the Jacobian:

A B

Jk= (14)

The flow term, Q, of the gravel equation [right-

∂R k ∂x

(23)

where k is the iteration step. Written in a compact matrix form stored as a seven-band matrix (suitable for band-matrix solvers; see Press et al., 1992, for example), the

T. den Bezemer et al. / Tectonophysics 320 (2000) 347–373

353

Jacobian matrix in Eq. (22) becomes: ∂R ∂R ∂R ∂R G,i−1 G,i−1 G,i−1 G,i−1 ∂h ∂F ∂h ∂F i−2 i−2 i−1 i−1 ∂R ∂R ∂R ∂R S,i−1 S,i−1 S,i−1 S,i−1 ∂F ∂h ∂F ∂h i−2 i−1 i−1 i ∂R ∂R ∂R ∂R G,i G,i G,i G,i ∂h ∂F ∂h ∂F i−1 i−1 i i ∂R ∂R ∂R ∂R S,i S,i S,i S,i ∂F ∂h ∂F ∂h i−1 i i i+1 ∂R ∂R ∂R ∂R G,i+1 G,i+1 G,i+1 G,i+1 ∂h ∂F ∂h ∂F i i i+1 i+1 ∂R ∂R ∂R ∂R S,i+1 S,i+1 S,i+1 S,i+1 ∂F ∂h ∂F ∂h i i i+1 i+1 i+2 We have coded the Newton iteration scheme by modifying a standard scheme by Press et al. (1992). This scheme is modified in four ways. First of all, a routine is supplied with an analytically computed Jacobian, rather than by using the finite difference Jacobian supplied by the standard routine. This makes the scheme considerably faster, and the Jacobian is easily obtained analytically (see below). Second, somewhat different stopping criteria are built in (see below). Third, a different matrix solver is used, which is especially designed for solving band matrices. To solve the matrix equations, we used the band-diagonal scheme in Press et al. (1992), which performs LU factorization and partial pivoting, when necessary. Fourth, we added a time-step chopping routine. When convergence was not achieved after a specified number of iterations, we divided the time step into 10 substeps and then solved the equation set for each of these substeps, which, in most cases, yielded convergence.

t 0 N N∂RS,i−1 N ∂hi−2 N 0 N N ∂h∂RS,i N i−1 N 0 N∂R v ∂hS,i+1

∂R G,i−1 ∂h i ∂R S,i−1 ∂F i ∂R G,i ∂h i+1 ∂R S,i ∂F i+1 ∂R G,i+1 ∂h i+2 ∂R S,i+1 ∂F i+2 method between

∂R G,i−1 ∂F i

u N 0 N N ∂R G,i N ∂F i+1 N. 0 N N ∂R G,i+1 N ∂F i+2 N 0 w

(24)

is to calculate an arithmetic average two nodal points.

F +F i±1 . F = i i±1/2 2

(25)

Another method to represent F and F i+1/2 i−1/2 is to realize that the composition of the incoming flux is dependent on surrounding high elevated columns, while the composition of the outgoing flux (to lower level columns) is determined only by the actual column (see Fig. 2). This approach is similar to the upwind differencing scheme in

5.6. Distance weighting of F The terms F and i+1/2 be evaluated properly. attempted to represent

F in Eq. (14) have to i−1/2 Two methods have been these terms. The easiest

Fig. 2. Visualizing upstream weighting of F. Column 3 receives sediment from (2). The gravel fraction that (3) receives depends on the gravel fraction of column (2). Similarly, column (4) is downstream to both (3) and (5), and hence, these two determine the gravel ratio of the flux calculated for (4).

354

T. den Bezemer et al. / Tectonophysics 320 (2000) 347–373

Press et al. (1992). The values for F and i+1/2 F are now governed by the height, h: i−1/2 F if h ≥h i i i±1 . F = i±1/2 F if h
G

(26)

The use of this upstream weighting of the fraction F makes the numerical solution of Eq. (2) somewhat more stable, and therefore, we have adapted this scheme in all the simulations in this paper.

5.7. Derivatives In the Jacobian, the following derivatives for the gravel equation have to be evaluated ∂M ∂Q ∂R G,i − G,i k=i−1, i, i+1 G,i = ∂h ∂h ∂h k k k ∂M ∂Q ∂R G,i − G,i k=i−1, i, i+1 G,i = ∂F ∂F ∂F k k k

∂Q K G,i = G ∂h Dx2 k ∂(hj+1 −hj+1 ) ∂(hj+1 −hj+1 ) i+1 i i−1 i +Fj+1 × Fj+1 i+1/2 i−1/2 ∂hj+1 ∂hj+1 k k k=i−1, i, i+1. (30)

A

B

Differentiated with respect to F: ∂Q

K G,i = G ∂Fj+1 Dx2 k ∂Fj+1 ∂Fj+1 i−1/2 (hj+1 −hj+1 )+ i+1/2 (hj+1 −hj+1 ) × i−1 i i+1 i ∂Fj+1 Fj+1 k k k=i−1, i, i+1. (31)

A

B

5.8. Initial guess for Newton iterations

(27)

(the same has to be done for the sand equation). The derivative of the accummulation term, M [Eq. (14)], depends on whether erosion or deposition has occurred. For erosion (hj+1
derivative differentiated with respect to h is:

(28)

Although we have been using a globally convergent scheme to find the root of the non-linear system of equations arising from the numerical discretization outlined above, it seems to be very sensitive to the initial guess for the root (e.g. Dennis and Schnabel, 1983; Press et al., 1992). We supply the multi-dimensional Newton solver with an initial guess, computed from an implicit Euler solution scheme for Eqs. (2) and (4), ignoring the non-linearities in F in the equation by taking Fj for Fj+1. This has shown to be very beneficial for reduction of the number of iterations needed for convergence (see Section 5.13).

With deposition (hj+1≥hj): i i

5.9. Scaling h and constraining F

∂M 3 3 ∂M 3 G = G = Fj+1 A+ hj+1 i ∂hj+1 2Dt ∂Fj+1 2Dt 2Dt i i i 1 4 hj + hj−1 . (29) − 2Dt i 2Dt i

The line-search procedure (Dennis and Schnabel, 1983; Press et al., 1992), ensuring global convergence in the applied Newton iteration scheme, assumes that typical values of the solution vector, x and R(x), are of order unity. To achieve this, we have rescaled the h-values, whose values are in the range of 0 to 1000 m, in the Newton scheme by dividing them by 1000 (note that A has to be scaled as well ). This rescaling improved the convergence properties. Another way of improving convergence was to constrain F between 0 and 1 from each iteration.

The derivative of the flow term, Q (15), depends on how the actual column is situated between two neighbouring columns, since upstream differencing of F is used. Therefore, we kept the F and i+1/2 F in the notation of the derivative. The partial i−1/2

355

T. den Bezemer et al. / Tectonophysics 320 (2000) 347–373

5.10. Stopping criteria Convergence is a matter of how the small time steps are compared to the rate of change within the system (Dennis and Schnabel, 1983). If the time step is too large, the routine needs too many iterations or may even never converge. We have been using two stopping criteria for convergence. The first is standard in the routine (Press et al., 1992) and checks for convergence of the solution vector, x. The second is an extension of the criterion supplied by the routine and checks whether

the Euclidean norm of the residual vector, R, is smaller than a predefined value and subsequently whether mass balance is achieved: 2N ∑ (R DxDt)2≤e. i i=1

(32)

5.11. Error decrease with decreasing time step The first analysis on numerical behaviour of the system compares the properties of the three implicit schemes in terms of accuracy (see Table 1).

Table 1 Analysis of error reduction with decreasing time step Dh

DF

fd time step (ka)

Number of time steps

h

F

Backward Euler (fully implicit) — erosion — K1=10 Ks t=100 t=50 t=25 t=12.5 t=6.25 t=3.125 t=1.50625 t=0.78125 t=0.390625 t=0.195312 t=0.097656

1 2 4 8 16 32 64 128 256 512 1024

102.85 92.58 86.63 83.51 81.93 81.15 80.75 80.56 80.47 80.41 80.37

0.9077 0.9067 0.9064 0.9064 0.9064 0.9064 0.9063 0.9063 0.9063 0.9063 0.9063

10.27 5.95 3.12 1.58 0.78 0.4 0.19 0.09 0.06 0.04 0

0.001 0.0003 0 0 0 1E−04 0 0 0 x 0

Crank–Nickelson — erosion — K1 t=100 t=50 t=25 t=12.5 t=6.25 t=3.125 t=1.50625 t=0.78125 t=0.390625 t=0.195312 t=0.097656

1 2 4 8 16 32 64 128 256 512 1024

– – – – 80.29 80.3 80.32 80.34 80.36 80.37 80.37

– – – – 0.8671 0.906 0.9063 0.9063 0.9063 0.9063 0.9063

– – – – −0.01 −0.02 −0.02 −0.02 −0.01 0 0

– – – – −0.0389 −0.0003 0 0 0 0 0

B3 (three time step method) — erosion — K1=10 Ks t=100 t=50 t=25 t=12.5 t=6.25 t=3.125 t=1.50625 t=0.78125 t=0.390625 t=0.195312 t=0.097656

1 2 4 8 16 32 64 128 256 512 1024

102.85 86.15 80.7 80.34 80.37 80.37 80.37 80.37 80.37 80.37 80.37

0.9077 0.9106 0.9065 0.9063 0.9063 0.9063 0.9063 0.9063 0.9063 0.9063 0.9063

16.7 5.45 0.36 −0.03 0 0 0 0 0 0 0

−0.0029 0.0039 0.0002 0.0002 0 0 0 0 0 0 0

356

T. den Bezemer et al. / Tectonophysics 320 (2000) 347–373

Fig. 3. Height-fraction profile used as input to test the numerical behaviour of the sedimentation model.

The decrease in Dh and DF is followed when the time step is halved. To analyse its behaviour, the system has been fed by a ‘difficult’ h–F profile (see Fig. 3). The implicit Euler scheme shows a linear decrease in Dh as expected (e.g. Smith, 1992). The Crank–Nicolson scheme has difficulties in obtaining a solution for larger time steps, whereas it does not converge quadratically for smaller time steps, as expected from a Crank–Nicolson scheme. The reason for this has not yet been investigated. The error, Dh, does decrease quadratically for the B3 method, as expected (Smith, 1992).

5.12. Comparison of solution schemes We have compared the four solution techniques (explicit Euler, implicit Euler, Crank–Nicolson, B3) also for consistency in results for both erosion and sedimentation. The results are shown in Table 2. For a time step of 0.1 ka, any differences are negligible. A simulation of 20 Ma duration for a scenario with changing subsidence rate (see Fig. 4) showed that also, visually, there are no detectable differences between the four solution techniques.

Table 2 Comparison of numerical solution techniques Time (ka)

Forward Euler h

Backward Euler F

h

F

Crank–Nickelson

B3

h

h

F

F

Time step=0.1 — erosion t=10 168.5 t=20 146.5 t=30 130.2 t=40 117.8 t=50 108.1 t=60 100.3 t=70 93.85 t=80 88.48 t=90 83.92 t=100 80.36

0.9051 0.9059 0.906 0.9062 0.9063 0.9063 0.9063 0.9063 0.9063 0.9063

168.6 146.5 130.3 117.9 108.2 100.3 93.92 88.55 83.97 80.36

0.9051 0.906 0.906 0.9062 0.9063 0.9063 0.9063 0.9064 0.9063 0.9063

168.6 146.5 130.3 117.8 108.1 100.3 93.88 88.52 83.95 80.37

0.9052 0.906 0.906 0.9062 0.9063 0.9063 0.9063 0.9064 0.9063 0.9063

168.6 146.4 130.1 117.7 108 100.2 93.83 88.45 83.9 80.37

0.9052 0.906 0.906 0.9062 0.9063 0.9063 0.9063 0.9064 0.9063 0.9063

Timestep=0.1 — sedimentation t=10 −135.9 t=20 −108.9 t=30 −93.15 t=40 −82.51 t=50 −74.72 t=60 −68.71 t=70 −63.91 t=80 −59.96 t=90 −56.64 t=100 −54.05

0.1625 0.2196 0.2626 0.2956 0.3217 0.3425 0.3594 0.3734 0.3851 0.3941

−136.2 −109.1 −93.3 −82.62 −74.81 −68.78 −64.01 −60.01 −56.69 −54.09

0.1621 0.2191 0.262 0.2951 0.3212 0.3421 0.3591 0.3731 0.3848 0.3939

−136.1 −109 −93.23 −82.56 −74.76 −68.74 −63.93 −59.99 −56.67 −54.07

0.1623 0.2193 0.2623 0.2954 0.3214 0.3422 0.3593 0.3733 0.3849 0.394

−135.7 −108.8 −93.1 −82.48 −74.69 −68.69 −63.89 −59.95 −56.64 −54.07

0.1629 0.2198 0.2626 0.2957 0.3216 0.3425 0.3594 0.3734 0.3851 0.394

T. den Bezemer et al. / Tectonophysics 320 (2000) 347–373

357

Fig. 4. Comparison of numerical solution schemes for a scenario involving a change in subsidence rate (see Section 6 for model set-up).

358

T. den Bezemer et al. / Tectonophysics 320 (2000) 347–373

equations is almost doubled for a solution without a calculated guess.

5.13. Convergence behaviour We have carried out a few tests to check and compare the convergence behaviour. We have compared the convergence behaviour for runs with and without a calculated initial guess (see Section 5.12). We have been feeding the routine the ‘difficult’ h–F profile again ( Fig. 3). The results are shown in Table 3. For large time steps, the number of iterations necessary to solve the set of

6. Perfect- and imperfect sorting for quasiequilibrium conditions Having discussed the numerical solution of the perfect-sorting and imperfect-sorting models above, predictions for a few basic scenarios involv-

Table 3 Analysis of convergence behaviour Iteration level

R1

R2

R3

R4

R5

R6

R

Convergence behaviour of Newton iterations — computed initial guess — t.d. time step=1 t=1 0 −0.068 −0.001 0 0 0 1 −2.50E−05 2.60E−04 0 0 7.20E−10 2 −7.20E−09 7.20E−08 0 0 1.20E−10 3 7.70E−14 −7.60E−13 0 0 1.30E−14 4 0 0 0 0 0

−1.50E−08 −7.20E−10 1.20E−09 −1.30E−13 0

1.20E+07 1.30E+05 232 0.17 0

Convergence behaviour of Newton iterations — computed initial guess — t.d. time step=0.1 t=1 0 −7.70E−02 −0.01 0 0 0 1 −2.20E−05 2.20E−04 0 0 0.00E+00 2 −5.60E−11 5.6E−10 0 0 0.00E+00 3 0.00E+00 0.00E+00 0 0 0.00E+00

0.00E+00 0.00E+00 0.00E+00 0.00E+00

1.40E+05 4.60E+03 9.3 0

Convergence behaviour of Newton iterations — initial guess=old values — t.d. time step=1 t=1 0 0 0.1 0 0 0 1 −7.70E−03 9.80E−02 0 0 −4.10E−07 2 −3.90E−03 4.90E−02 0 0 −2.10E−07 3 −2.00E−05 2.00E−04 0 0 −8.30E−13 4 8.60E−10 −8.60E−09 0 0 6.80E−15 5 −3.40E−16 3.40E−15 0 0 0 6 0 0 0 0 0 7 0 0 0 0 0 8 0 0 0 0 0 9 0 0 0 0 0 10 0 0 0 0 0

0 4.10E−07 2.10E−07 −8.30E−12 −6.80E−14 0 0 0 0 0 0

1.60E+07 1.50E+06 1.40E+06 6.30E+05 4.50E+05 4.40E+05 2.60E+05 2.50E+05 6.70E+04 5.20E+04 0

Convergence behaviour of Newton iterations — computed initial guess — t.d. time step=1 t=10 0 2.60E−04 1.70E−05 0 −1.70E−08 2.80E−05 1 6.30E−08 −1.00E−07 2.30E−11 −2.40E−10 −1.70E−09 2 −8.20E−15 −4.40E−13 0 0 1.50E−16 3 0.00E+00 0.00E+00 0 0 00E+00

−2.80E−05 −4.50E−11 −3.10E−15 00E+00

6.50E+04 5.10E+01 0.015 0

Convergence behaviour of Newton iterations — initial guess=old values — t.d. time step=1 t=10 0 7.60E−03 3.503−05 1.90E−07 −2.00E−08 8.10E−10 1 4.30E−06 5.80E−06 9.10E−10 9.30E−09 2.50E−10 2 8.50E−11 −4.90E−13 −7.60E−14 7.50E−14 6.00E−15 3 00E+00 00E+00 0 0 00E+00

2.30E−05 2.10E−10 00E+00 00E+00

1.40E+05 1.03E+03 2.7 0

SUM

T. den Bezemer et al. / Tectonophysics 320 (2000) 347–373

ing periodic changes of variables governing the diffusion equations will be compared. Based on analogues from heat-conduction theory, Angevine et al. (1990) presented a non-dimensional analysis of the basic topographic diffusion equation. Adopting the characteristic scales, h , t , and L 0 0 for h, t and x, respectively and denoting nondimensional variables with an asterisk, Eq. (1) transforms to: h ∂h1 kh ∂2h1 0 = 0 t ∂t1 L2 ∂x21 0 and upon rearranging: L2 ∂h1

=

(33)

∂2h1

. (34) kt ∂t1 ∂x21 0 The dimensionless number, L2/k, denotes the natural timescale for the basin, and if the timescale of variations imposed on the basin is long compared to this natural scale, the basin remains in quasiequilibrium as it responds to these changes (Angevine et al., 1990). Slow variations and fast variations are then defined as variations on timescales greater than, and less than, L2/k, respectively. Calculating the natural timescale for the dual lithology diffusion equations is not as simple as it is for the basic topographic diffusion equation. The problem is in calculating the effective diffusion coefficient [see height equation, Eq. (3)]. Consequently, it is also very difficult to match the diffusion coefficient of the basic topographic diffusion equation with an effective diffusion equation in the dual lithology system. In fact, the dual lithology system produces a somewhat different depositional profile, due to the advection term present in the height equation. Writing out Eq. (3):

A B B

359

system. This makes any precise comparison rather difficult. Being interested in first-order comparisons, we chose a very straightforward approach and simply chose diffusion coefficients, K and G K , in such a way that their arithmetic mean was S equal to k in the simple topographic diffusion equation (see below). We will start by presenting model results for quasi-equilibrium variations, i.e. variations on timescale long compared with the natural timescale of the basin. The basin length is 20 km, and the diffusivity is 4×102 m2/a in the perfect-sorting model. The response of the basin L2/k then equals 106 years. The diffusivities in the imperfect-sorting model are 6×102 m2/a and 2×102 m2/a for K G and K , respectively. Sinusoidal variation in subsiS dence rate, gravel fraction and diffusivity (which is mainly controlled by water supply; see Appendix A) occurs on a timescale of 107 years, which is 10 times the natural response time. As in the simulations by Angevine et al. (1990), the basin has a simple asymmetric distribution of subsidence increasing linearly to the left (see Fig. 5). For both sedimentation models, constant height boundary conditions have been used. Note that by taking fixed height boundary conditions (as we have done throughout this paper), the sediment input is controlled by the subsidence rate and diffusivity. A fixed height boundary is a natural situation for locally controlled fault-related sedimentation. Coupling of sediment input with subsidence and diffusivity can be avoided, however, by taking other, fixed flux-type boundary conditions (Angevine et al., 1990; Rivenaes, 1992).

∂2h =((K −K )F+K ) G S S ∂t ∂x2

∂h

A

∂F + (K −K ) G S ∂x

∂h

∂x

.

(35)

This is a diffusion–advection equation (for its behaviour, see Kooi and Beaumont, 1994). The diffusion term has a tendency to smooth a topographic profile, whereas the advection term propagates a part of a topographic profile through the

Fig. 5. Cross-section showing simulated basin conditions. Black arrows are subsidence arrows and grey arrows are sedimentation arrows (after Angevine et al., 1990).

360

T. den Bezemer et al. / Tectonophysics 320 (2000) 347–373

Fig. 6. Cross-section showing a hypothetical basin under the influence of a slow variation in subsidence rate. Figure conventions: the upper panel shows the predicted grain-size distribution using a perfect-sorting model, and the lower panel shows the results from an imperfect-sorting model. The used parameters are shown below.

Composition of sediment input to the basin for the perfect-sorting model has been done using a fixed gravel flux. In all simulations with the imper-

fect-sorting model, the right-hand side fraction is adjusted in such a way that it delivers a constant fraction of 0.5 to the basin.

T. den Bezemer et al. / Tectonophysics 320 (2000) 347–373

This has been done by setting: K S . (36) K +K S G The simulation time for one of the equilibrium models is typically a few hours on a Sun UltraSparc.

F=

6.1. Slow variation in subsidence The first model run compares predictions of the perfect- and imperfect-sorting models for a slowly varying subsidence. Because we are using fixed height boundary conditions, a varying subsidence means a varying sediment boundary flux as well. Sedimentation rates increase as the subsidence increases (see Fig. 6). The perfect-sorting model (upper panel in Fig. 6) predicts progradation of gravel during periods of reduced subsidence when the proximal accommodation space for gravel is small. The imperfect-sorting model also predicts coarse progradation during times of reduced subsidence on the left-hand side of the basin. The righthand side of the basin, however, shows exactly the reverse; during phases of increased subsidence, progradation occurs. This behaviour is related to the asymmetry of subsidence. With increasing subsidence rate, the basin axis tends to shift to the left (Fig. 7). This causes coarse fractions coming from the left flank of the model to become trapped on the left-hand side of the basin. At the same time, sediments from the right-hand side of the basin are transported more towards the left in the direction of the shifting basin axis. For decreasing subsidence, the opposite occurs.

Fig. 7. Cross-section showing changes in basin shape under influence of changing diffusivity or subsidence rate.

361

6.2. Slow variation in diffusivities A slow variation in diffusivity has no visible effect on either the sedimentation rate or the gravel front position in the perfect-sorting model ( Fig. 8, upper panel ). This model behaviour is similar to the equivalent simulation by Angevine et al. (1990) and is explained by the fact that the diffusivity is not an important element in the long-term balance of Eq. (34). In contrast with the perfect-sorting model, the imperfect-sorting model shows (Fig. 8, lower panel ) pronounced shifting facies (boundary fractions are adjusted to deliver a constant gravel input to the basin while changing the diffusivities of gravel and sand). Similar to the perfect-sorting model, there are no visible changes in sedimentation rate. Shifting of facies is caused by a changing basin shape. As in the previous simulation (varying subsidence), the grain-size pattern on the righthand side of the basin shifts out of phase with that of the left-hand side of the basin. This is due to the same competition-of-slope effect discussed in the Section 6.1. During periods of high diffusivities, the left-hand slope is dominant, whereas during periods of low diffusivities, the right-hand slope is dominant. Note that the sorting efficiency does not change as a result of varying both K and K with G S amplitudes of 70% of their mean values. Although the absolute value of the difference between K G and K changes ( Fig. 9), their ratio — which is S important — does not. The next simulation deals with the situation in which this ratio does change.

6.3. Slow variation in diffusivity gravel only The lower panel of Fig. 10 shows the model behaviour when only the diffusivity of gravel is slowly varied. Initially, the diffusivity ratio between the gravel and sand decreases, which results in a very broad zone between the coarsest (purple colours) and finest (blue–green colours) sediments. When the diffusivity ratio starts to rise again, a very narrow zone between the coarsest and finest sediments is developed as a result of an increased sorting efficiency. Varying this K on long timesG cales does not affect sedimentation rate much, but,

362

T. den Bezemer et al. / Tectonophysics 320 (2000) 347–373

Fig. 8. Cross-section showing a hypothetical basin under the influence of a varying diffusivity. Figure convention as in Fig. 6.

as in the previous experiment, it does influence the position of the basin axis. During periods of high gravel diffusivity, the total effective diffusivity is

raised, and the basin axis shifts to the right, causing coarse facies on the left to prograde and coarse facies on the right to retrograde.

T. den Bezemer et al. / Tectonophysics 320 (2000) 347–373

363

Fig. 9. Changes in the difference in diffusivity between gravel and sand.

6.4. Slow variation in gravel supply (fraction supply) The last slow variation considered is the variation in gravel supply (fraction). This is probably the most straightforward simulation to interpret in both the perfect- and imperfect-sorting model. A change in gravel fraction supplied to the basin is likely to occur when the source rock lithology changes, due to faulting or unroofing or when climate changes and reduces the generation of fines by weathering. Both models show shifting of facies without any changes in accumulation rate (see Fig. 11). Note, also, that although the fraction delivered to the basin by the right-hand boundary is kept constant, it shifts exactly in phase with the facies on the left-hand side (Fig. 11, lower panel ). This is again a competition of slope effects ( Eq. 26 and Fig. 7). During phases of low fraction input (few gravel ), the sedimentation rate is governed by higher diffusing sands, so the left-hand slope

will be relatively dominant compared to the phases of high gravel input.

7. Perfect- and imperfect sorting for nonequilibrium conditions The simulations in this section were made under the same conditions as in Section 6, except that variations now occur on a timescale of a factor 10 less than the equilibrium time. A model run takes typically less than 30 min on a Sun UltraSparc. In contrast to the previous simulations, the righthand side of the basin is not discussed, mainly because it is barely visible in the figures. 7.1. Fast variation in subsidence rate The fast variation in subsidence shows gravel progradation during phases of reduced subsidence

364

T. den Bezemer et al. / Tectonophysics 320 (2000) 347–373

Fig. 10. Hypothetical basin when the input of gravel the fraction is varied. Figure convention as in Fig. 6.

T. den Bezemer et al. / Tectonophysics 320 (2000) 347–373

365

Fig. 11. Resulting pattern when only one diffusivity is varied. The upper panel shows the resulting pattern from a fast variation in diffusivity, and the lower panel shows the result of a slow variation in diffusivity. Figure convention as in Fig. 6.

366

T. den Bezemer et al. / Tectonophysics 320 (2000) 347–373

Fig. 12. Cross-section showing a hypothetical basin formed under the influence of a fast variation in subsidence rate. Figure convention as in Fig. 6.

T. den Bezemer et al. / Tectonophysics 320 (2000) 347–373

(see Fig. 12) in much the same way as it does for a slow variation in subsidence. Both models predict an outward building of coarse facies during phases of reduced subsidence. 7.2. Fast variation in diffusivities When diffusivity is varied on a short timescale (Fig. 13), both models predict very different grainsize distributions. The perfect-sorting model predicts coarse outbuilding when the diffusivity is low. In contrast, the imperfect-sorting predicts outbuilding when the diffusivity is high. This simulation clearly shows the difference in the modelling approach of both sedimentation models. The perfect-sorting model uses a constant absolute amount of gravel input on the left-hand side. During periods of decreasing k (transport efficiency), the mean sedimentation rate of the section decreases. Also, on the left-hand side of the model, sedimentation rates decrease. To accommodate the same volume of gravel per unit time, outbuilding/progradation has to occur. For increasing k, the opposite occurs. Shifts in grain-size boundaries in the imperfectsorting model are not caused by temporal changes in the fraction flux (the fraction, F, on the lefthand side is constant; only the magnitude of the flux is variable). In this model, facies shifts are induced by lateral migration of the basin axis. During decreasing k, the basin axis moves to the left, and facies are compressed towards the lefthand side of the model. For an increasing k, the opposite occurs. 7.3. Fast variation in diffusivity gravel only The upper panel of Fig. 10 shows the model behaviour when only one diffusivity, the diffusivity of gravel, is varied. In contrast with the slow equivalent of this simulation ( lower panel in Fig. 10), a fast variation in the diffusivity does change the accumulation rate, which results in erosional surfaces during periods of low diffusivity. Like the slow simulation, during periods of low gravel diffusivity, the horizontal distance between the coarse and fine fraction sediment is reduced.

367

7.4. Fast variation in gravel supply (fraction supply) A fast variation in gravel supply gives rather similar predictions for perfect- and imperfect sorting ( Fig. 14). As in the quasi-equilibrium situation, the perfect-sorting model predicts outbuilding of gravel without any changes in accumulation rate. The imperfect-sorting predicts clean peaks of coarse facies entering the basin, but unlike the perfect-sorting model, it predicts small changes in the accumulation in the centre of the basin. This is due to the fact that, during phases of high coarse input, the basin is dominated by a high fraction in the surface layer, which lowers the effective diffusivity (coarse facies has lower diffusivity) and hence decreases transport to the basin centre.

8. Discussion and conclusions During the last 30 years, many techniques have been developed to predict grain-size distributions in the stratigraphic record. Many of these models use some kind of equilibrium assumption, either for the depositional profile or for the relation between a certain grain size and depositional slope. In this paper, two models were analysed, a perfectsorting model and an imperfect-sorting model, which do not make these assumptions. The numerical solution of one these models, the imperfect-sorting model, is not straightforward. We have come up with a practical and robust scheme, the three-time-step method. The convergence properties for the Newton scheme necessary to solve the non-linear partial differential equations are considerably improved by calculating an initial guess, using an implicit Euler scheme using F(raction)-values of the previous time step. The model behaviour of the perfect- and imperfect-sorting model is rather similar for slowly (relative to the natural timescale of the basin) varying rates of subsidence and gravel input. The grain-size distribution, as predicted by the perfect- and imperfect-sorting models, however, is controlled by very different parameters. Gravel progradation in the perfect-sorting model, in the absence of gravel input variations, is controlled by

368

T. den Bezemer et al. / Tectonophysics 320 (2000) 347–373

Fig. 13. Cross-section showing grain-size distribution for a basin under the influence of a fast variation in diffusivity. Figure convention as in Fig. 6.

T. den Bezemer et al. / Tectonophysics 320 (2000) 347–373

Fig. 14. Hypothetical basin as a result of a fast variation in the input of the gravel fraction. Figure convention as in Fig. 6.

369

370

T. den Bezemer et al. / Tectonophysics 320 (2000) 347–373

all processes that decrease the accommodation space at the side of gravel input and, of course, by the amount of gravel input. Gravel progradation occurs when the subsidence rate decreases, the transport efficiency (k) decreases, or the gravel input increases. Shifting of facies in the imperfect-sorting model is more complicated. Shifts are controlled by the basin geometry, sorting efficiency (difference between diffusivities of both grain-size classes) and input fraction. Progradation of coarse facies occurs when the basin axis shifts away from that side of the basin and when the sorting efficiency is greatest. A basin axis shift may occur as a result of changes in subsidence rate or diffusivity (transport efficiency, water supply). The variation in diffusivity yielded the most complex model response, and the two models responded very differently to it. This is because the variation in diffusivity in the imperfect-sorting model affects the position of the basin axis in both equilibrium and non-equilibrium situations, causing a shifting of facies not observed in the perfectsorting model. This observation suggests that gravel progradation without any changes in accumulation rate can be caused by both changes in gravel supply and long-term changes in water supply (diffusivity) (see Figs. 10 and 8). Changes in the width of the transition zone between coarse and fine facies are caused by changes in the sorting efficiency (Fig. 10) and may lead to proximal coarsening, contemporaneous with distal fining. The perfect-sorting model suggests that gravel progradation as a result of changes in subsidence rate may be distinguished from gravel progradation resulting from changes in supply, using accumulation rates ( Figs. 12 and 14). The imperfectsorting model, however, shows that changes in accumulation rate may also occur as a result of (fast) changes in supply (Fig. 14), making this criterion not completely reliable. When judging the implications for interpretation of the geological record, one has to decide which model describes sediment transport in alluvial basins most realistically. Although perfect sorting is believed to occur in alluvial fans as a result of protrusion and hiding effects (Paola, 1988), imperfect sorting is clearly much more

realistic (Rivenaes, 1992). If imperfect sorting is the more realistic selective deposition mechanism, modelling results suggest that interpretation of progradation of coarse facies has to be done cautiously. Progradation may occur as a result of (1) increased coarse input, (2) lower subsidence rates, (3) higher subsidence rates (right-hand side basin, Fig. 6), (4) higher transport efficiency (water supply), and (5) lower transport efficiency (see Fig. 8, right-hand side basin). Differentiation between these controls can only be done using timelines and basin shapes.

Appendix A: Derivations of the sedimentation models This section largely follows Angevine et al. (1990). The simplest basis on which a dynamic model of alluvial deposition can be constructed is the diffusion equation. The diffusion equation has been used as a metaphor for sedimentation and erosion in a number of different sedimentary environments ( Kenyon and Turcotte, 1985; Begin, 1987; Angevine et al., 1990; Kooi and Beaumont, 1994). The derivation given below is based on the derivations by Begin (1987) and Angevine et al. (1990). Using conservation of mass and momentum for both water and sediment, the linear diffusion equation will be derived. Consider a two-dimensional, vertical slice of a sedimentary basin. Two-dimensional in this case means that the rates of subsidence and sediment supply do not vary perpendicular to the direction of transport. The total water discharge, q (m2/s), may be expressed as: q=bgu

(A1)

where q is the total water discharge per unit width of basin, b (−) the total width of channels per unit width of basin, g (m) the mean channel depth and u (m/s) the mean velocity. For spatial scales that are long compared with those of the topography on the river system (e.g. bars and channel bends), the conservation of

371

T. den Bezemer et al. / Tectonophysics 320 (2000) 347–373

momentum may be written as: t=−grg

∂h

(A2)

∂x

where t (N/m2) is the kinematic shear stress, g (m/s2) the gravitational acceleration, r (kg/m3) the density of the water, h (m) the elevation of the sediment relative to a non-subsiding horizontal datum and x (m) the downstream distance. Conservation of mass for sediment may be expressed as: ∂h

1 ∂( b q ) s (A3) ∂t C ∂x 0 with s (m/s) subsidence rate, t (s) the time, C (−) 0 the volume concentration of sediment in the bed, and q (m2/s) the rate of sediment transport per s unit width of channel. The angle brackets around q indicate that this is a long-term average. s Only a few semi-empirical relations are available for mass transport (Angevine et al., 1990).

s=

(A4)

t (N/m2) is the critical shear stress needed to set c the sediment in motion, and s is the sediment specific gravity (−). This relationship applies for short timescales. In order to combine this equation with the long-term equations, one has to assume that the flow can be characterized roughly by a set of representative channel-forming conditions that occur intermittently (Angevine et al., 1990). This leads to:

q =Iq . (A5) s s where I (−) is an intermittency factor. Combining Eqs. (A3) and (A4) yields: ∂(b(t−t )3/2) c . (A6) ∂t C g(s−1) ∂x 0 The shear stress can be easily related to velocity: ∂h

(A8) c where e is a constant. Introducing Eq. (A8) into Eq. (A6) gives:

t=(1+e)t

=

8(t−t )3/2 c q= . s g(s−1)

s=

Parker (1978) found a relationship that related actual shear stress to critical shear stress. He demonstrated that there was a self-regulating mechanism for the stress in the channel centre. If stress in the channel centre rises above a certain value, the banks begin to erode, thereby widening the channel and reducing the stress in the centre. Stress is maintained as:

=

−8I

t=c u2 (A7) f where c (kg/m) is the drag coefficient and is f usually taken to be constant. For gravel-bed rivers with non-cohesive banks,

∂h

−8I =

A

e (1+e)

B

3/2 ∂( bt3/2)

. C g(s−1) ∂x 0 Using Eqs. (A1), (A2) and (A7) yields:

s+

s+

∂t

∂h ∂t

=k

∂2h ∂x2

(A9)

(A10)

where k, the diffusivity is given by: −8 qA Ec

f. (A11) C (s−1) 0 Paola (1988) makes use of this relationship to find constraints for the value of k. I have not done this recognizing the enormous variability in reported values ( Kooi and Beaumont, 1994). The diffusion equation [Eq. (A10)] is well known from heat conduction problems (e.g. Carslaw and Jaeger, 1973), and its solution is only briefly discussed in chapter 2.

k=

Derivation of the dual lithology equation This part of the appendix largely follows Rivenaes (1992). Conservation of mass requires that sediment that accumulates in a certain control volume, during a short time, is equal to the difference between the mass that enters and leaves the control volume. The adjacent figure (A1) shows gravel fluxes (q , q ) entering and leaving this volume. G,in G,out Inside the volume, the gravel constitutes a bulk fraction, F, while the rest of the volume (1−F ) is occupied by sand. Within the gravel, C is the G

372

T. den Bezemer et al. / Tectonophysics 320 (2000) 347–373

and gravel solid fraction, C , in the surface layer G will be almost invariant. Dividing by DtDx on each side and letting Dt and Dx tend to zero results in: AC

∂h ∂q ∂F +C F =− G G ∂t G ∂t ∂x

(A14)

The average term C F should be treated somewhat G differently, depending on whether erosion of deposition occurs at a given position. With erosion, the C and F of the removed layers are known G and can be calculated using:

Fig. A1. Gravel entering and leaving the control volume (after Rivenaes, 1992).

solid fraction, that is (1−C ) is the porosity. In G order to simulate armouring, a surface transport layer with thickness, A, and bulk fraction, F, is used. This layer represents the conduit for lateral migration of sediments, containing all processes that transport and sort sediment (channel transport, creep, slumps and slides). For simplicity, it is assumed to be constant over the whole profile and independent of time. The thickness of the surface layer cannot be zero, so there should always be soil available for erosion. The mass balance for gravel is: ([C Fr ]t+Dt−[C Fr ]t)ADx G G G G +C Fr Dx(ht+Dt−ht)=[q r Dt]x G G G G −[q r Dt]x+Dx. (A12) G G Superscripts t and t+Dt represent initial time and new time, respectively. The first term on the lefthand side represents the net mass loss/gain within the boundary later, and the second term represents the mass loss/gain beneath the boundary layer. The right-hand side represents net mass loss due to transport. A precise expression for the averages in the third term is given by: 1

[C Fr ]= G G ht+Dt−ht

P

ht+Dt−A ht−A

(C Fr ) dy. G G

1 [C F ]= G ht+Dt−ht

P

ht+Dt−A

(C F ) dy. G

(A15)

ht−A If deposition occurs, the value of C F is unknown G and must be a part of the solution; F is equal to the F in the surface layer, and C is set to the G initial value. [C F ]=Ct+Dt Ft+Dt. (A16) G G The gravel volume rate, q , has to be replaced by G a sediment transport relation. This relation assumes that the flux, q, is a function of the slope. q=−K

∂h ∂x

.

(A17)

The gravel volume rate, q , must be a fraction of G the total flux of sediments. This fraction is not known, but intuitively, two requirements must be fulfilled: the fraction is zero when no gravel is available in the boundary layer, and conversely, if only gravel is available, the gravel flux equals the total flux. If a linear relationship is assumed and one realizes that K is lithology-dependent, the gravel equation is expressed by: AC

∂F

A

B

∂h ∂h ∂ FK . +C F = G ∂t G ∂t G ∂x ∂x

(A18)

The sand equation can be derived similarly:

A

B

∂h ∂h ∂F ∂ FK . −AC +C (1−F ) = S G S ∂t ∂t ∂x ∂x (A19)

(A13) It is reasonable to suppose that the grain density

In this dual-lithology system, the set of non-linear coupled partial differential equations has to be

T. den Bezemer et al. / Tectonophysics 320 (2000) 347–373

solved for the two unknowns in the system, F and h. The solution is shown extensively in Appendix A.

References Angevine, C.L., Heller, P.L., Paola, C., 1990. Quantitative sedimentary basin modelling. AAPG Continuing Education Course Note Series Vol. 32. Begin, Z.B., 1987. Application of a diffusion–erosion model to alluvial channels which degrade due to baselevel lowering. Earth Surf. Proc. Landforms 13, 487–500. Brierley, G.J., Hickin, E.J., 1985. The downstream gradation of particle sizes in the Squamish River, British Columbia. Earth Surface Processes Landforms 10, 597–606. Carslaw, H.S., Jaeger, J.C., 1973. Conduction of Heat in Solids. second ed., Oxford University Press. Cross, T., Harbaugh, J., 1990. Quantitative dynamic stratigraphy: a workshop a philosophy a methodology. In: Cross, T.A. (Ed.), Quantitative Dynamic Stratigraphy, 625. Dennis, J.E., Schnabel, R.B., 1983. Numerical Methods for Unconstrained Optimization and Non-linear Equations. Prentice Hall, Englewood Cliffs, NJ, 378 pp. Ferguson, R., Hoey, T., Wathen, S., Werritty, A., 1996. Field evidence for rapid downstream fining of river gravels through selective transport. Geology 24, 179–182. Flemings, P.B., Grotzinger, J.P., 1996. STRATA: freeware for analyzing classic stratigraphic problems. GSA Today 6 (12), 1–7. Hardy, S., Waltham, D., 1992. Computer modeling of tectonics, eustacy, and sedimentation using the Macintosh. Geobyte 7 (6), 42–52. Jervey, M.T., 1988. Quantitative geologic modeling of siliclastic rock sequences and their seismic expression, Wilgus, C.K., Hastings, B.S.Posamentier, H., Wagoner, J.V., Ross, C.A., Kendall, C.G.St.C. ( Eds.), Sealevel Changes — An Integrated Approach. SEPM Special Publication 42, 47–69.

373

Kenyon, P.M., Turcotte, D.L., 1985. Morphology of a delta prograding by bulk sediment transport. Geological Society of America Bulletin 96, 1457–1465. Kooi, H., Beaumont, C., 1994. Escarpment evolution on highelevation rifted margins: Insights derived from surface processes model that combines diffusion, advection and reaction. Journal of Geophysical Research 99, 12 191–12 209. Krumbein, W.C., 1941. The effects of abrasion on the size, shape and roundness of rock fragments. Journal of Geology 49, 482–520. Lawrence, D.T., Doyle, M., Aigner, T., 1990. Stratigraphic simulation of sedimentary basins: concepts and calibration. AAPG Bulletin 74 (3), 273–295. Paola, C., 1988. Subsidence and gravel transport in alluvial basins. In: Kleinspehn, K.L., Paola, C. ( Eds.), New Perspectives in Basin Analysis. Springer, Berlin, p. 453. Parker, G., 1978. Self-formed straight rivers with equilibrium banks and mobile bed, Part 2, the gravel river. Journal of Fluid Mechanics 89, 127–146. Pitman, W.C., 1978. Relationship between eustacy and stratigraphic sequences of passive margins. Geological Society of America Bulletin 89, 1389–1403. Press, W.H., Flannery, B.P., Teukolsky, S.A., Vetterling, W.T., 1992. Numerical Recipes in C. Cambridge University Press, Cambridge, 994 pp. Rivenaes, J.C., 1992. Application of a dual-lithology, depthdependent diffusion equation in stratigraphic simulation. Basin Research 4, 133–146. Shaw, J., Kellerhals, R., 1982. The composition of recent alluvial gravels in Alberta River beds: Alberta. Research Council Bulletin 41, 151. Shumm, S.A., Stevens, M.A., 1973. Abrasion in place: a mechanism for rounding and size reduction of coarse sediments in rivers. Geology 1, 37–40. Smith, G.D., 1992. Numerical Solution of Partial Differentiation Equations. Oxford University Press, Oxford, 337 pp. Strang, G., 1988. Linear Algebra and its Application. Harcourt Brace Jovanovich, 505 pp. Tetzlaff, D.M., Harbaugh, J.W., 1989. Simulating Clastic Sedimentation. Van Nostrand Reinhold, New York, 202 pp.