7.3.1 Dynamics of a Single Soil Layer

We will first describe how mARM3D simulates each layer, and then we will discuss how the layers interact vertically. We start with the thin surface layer. This layer is the layer that interacts with the water flowing over the surface. In the most general case if the sediment transport capacity of the flow is more than the amount of sediment being carried by the flow, then the flow will erode material from the surface, while if the transport capacity is less than the sediment in the flow, then it will deposit sediment. If erosion occurs, then there is preferential entrainment of the finest fractions from the surface layer into the flow. To maintain the mass of the surface layer when there is erosion, material is transferred from the layer directly below to exactly balance the material being removed by erosion. If deposition occurs, then there is preferential deposition of the coarsest fractions (the coarsest particles settle out fastest) from the flow into the surface layer. To maintain the mass balance in the surface layer during deposition, material from the surface layer is transferred into the layer directly below to balance that material being deposited.

For the layers below the surface layer the mass balance is maintained at every time step. Thus if material is transferred from the layer to be put into another higher layer, then material is entrained from the layer below to balance the mass removed. Likewise if material is transferred into the layer from the layer above, then material from that layer is pushed down into the layer below. In the simplest case, if material is removed from the top layer by erosion, then material from the next layer below is transferred up to balance the mass lost to erosion. The next layer below that then has material transferred upwards to balance the material lost from the layer above, and this processes cascades down the profile to the bottom of the soil profile so that the soil-bedrock interface moves closer to the surface and the soil thins.

Within each layer weathering can occur and that weathering changes the grading distribution of that layer. Each layer weathers independently of the other layers. Mass conservative (i.e. physical weathering) and non-mass conservative (i.e. dissolution) weathering can be modelled, but to date full chemical weathering where dissolution is followed by precipitation of secondary minerals has not been modelled, because mARM does not have a coupled biogeochemical model to simulate in-profile geochemistry (Chapter 8). The weathering rate is typically depth dependent and can also be a function of particle grading and time. The model does not have a coupled groundwater model and so cannot model the effects of the interaction between soil moisture, soil grading and weathering rates, though a known, specified soil moisture spatial pattern (and thus a specified spatial weathering pattern) can be input.

The bottom layer of the soil profile has the grading of the underlying rock, rock being defined as material that is 100% the coarsest size fraction in the particle size distribution. That layer becomes soil as the weathering function breaks down the material into smaller size fractions. Any layers below cannot be soil until a specified proportion of the rock fraction in the layer above is broken down. Thus the bedrock-soil interface arises naturally as a result of the weathering process and is not explicitly modelled by the soil production function from Chapter 6.

Having outlined the conceptual approach used in mARM3D the mathematical details follow. The approach draws heavily from the state-space literature, and some of the terminology of this literature will be used where appropriate.

The soil grading at any given time and in any given layer is represented by a vector called the state vector $\underset{¯}{g}:$

$\underset{¯}{g}=\left(\begin{array}{c}{g}_1\\ g_2\\ ⋮\\ g_{m-1}\\ g_m\end{array}\right)$(7.2)

and the entries in the state vector, g_i,$g_i,$ are the mass of sediment in each grading size range i where there are m size fractions for the grading. Cohen et al. (Reference Cohen, Willgoose and Hancock2010) used the proportion by mass of the layer in each size range as the state, but subsequent experience has found that using the actual mass in the layer in each size range makes it easier to apply mass conservation principles directly to the construction of the transition matrices (see below). The transition from the grading at any given time to the grading at the next time step is described by a matrix equation. It describes both how the grading changes with time and how the grading at one time, t₁, is related to the grading at some time in the future t₂:

$\underset{¯}{g}\left(t_2\right)=R\underset{¯}{g}\left(t_1\right)$(7.3)

where in the state-space literature the matrix R is called the transition matrix and is typically a function of the size of the timestep (t₂ − t₁).$\left(t_2-t_1\right).$ Note that any set of coupled differential equations can also be expressed with Equation (7.3). The advantage of the matrix formulation is that we can explicitly formulate all the processes that change the grading with the same methodology.

To represent the soil profile, we write one of these equations for each layer so that we can write Equation (7.3) for each layer down the profile, where R_j is different in each layer reflecting how weathering processes change down the profile. For each layer j we then have

${\underset{¯}{g}}_j\left(t_2\right)=R_j{\underset{¯}{g}}_j\left(t_1\right)$(7.4)

The discussion in the remainder of this section will consider only what happens in a single layer, so we will drop the j subscript for the moment.

We simulate each physical process as a multiplicative change to the state. If we have a single process, let’s call it A, with a corresponding transition matrix R, then the evolution of the grading vector over one timestep is (using subscripts for the timestep)

${\underset{¯}{g}}_{t+1}=R{\underset{¯}{g}}_t$(7.5)

which is simply Equation (7.3) where one timestep is (t₂ − t₁).$\left(t_2-t_1\right).$ To demonstrate how the method works for more than one process, consider the case where there are two independent physical processes, called A and B with corresponding transition matrices in Equation (7.5) of R and S; then the combination of these processes on the soil grading is

${\underset{¯}{g}}_{t+1}=S\left(R{\underset{¯}{g}}_t\right)=\mathrm{SR}{\underset{¯}{g}}_t$(7.6)

where the order of the matrix multiplication implies that process A acts on the grading first and B operates second on the result of process A. This can be generalized to any number of processes, and is the key to the simplicity and generality of the approach. This splitting of processes into independent operations is ‘operator splitting’ (Section 2.2.2), and this splitting considerably simplifies the task of constructing the matrices for the combined effect of simultaneous physical processes.

Equations (7.3) to (7.6) follow the traditional presentation of transition matrices. In the discussion that follows it is more convenient to work with the marginal transition matrix than the actual transition matrix. In this case

$R=I+A$(7.7)

where I is the identity matrix (where $\underset{¯}{g}=I\underset{¯}{g}$) and A is the marginal transition matrix for process A. The change in grading over a timestep is proportional to matrix A. The advantage of formulating the matrices in this form is that the rate of the process represented by A (e.g. the weathering rate) can be changed simply by multiplying the marginal matrix by a scaling factor. For example, doubling the process rate is achieved by applying Equation (7.5) twice (the same as doing two timesteps), and if the timestep is small, then this is the same as multiplying the marginal transition matrix by 2 so that

$\begin{array}{c}\text{process}\text{A}\left(\text{nominal rate}\right)=R=I+A\\ \text{process}\text{A}\left(\text{twice nominal rate}\right)=R^2={\left(I+A\right)}^2\cong I+2A\end{array}$(7.8)

where the A²$A^2$ term is dropped since it is small for a small timestep. In the discussion that follows, the marginal transition matrix is used unless otherwise noted. In marginal transitional matrix form, Equation (7.6) is

${\underset{¯}{g}}_{t+1}=\mathrm{SR}{\underset{¯}{g}}_t=\left(I+B\right)\left(I+A\right){\underset{¯}{g}}_t$(7.9)

where A and B are the marginal transition matrices for processes A and B, respectively, and correspond to the transition matrices R and S in Equation (7.6). Representing weathering and other processes within a layer is then a matter of constructing the marginal transition matrices and repeatedly applying the equations above.

Before we discuss how soil layers interact, a simple example will be useful to explain how the details of this process work. Consider an example where there are only three grading size ranges (for convenience let’s call them i = 1,$i=1,$ 2, 3, ‘small’, ‘medium’ and ‘large’ grain sizes), and we will construct a marginal transition matrix for a weathering process that breaks large particles into medium particles, medium particles into small particles, and leaves the small particles unchanged. If at any one timestep 2% of the mass of the large particles are weathered into medium sized particles and 1% of the medium into fine, then the marginal transitional matrix is

$A=\left(\begin{array}{ccc}0& 0.01& 0\\ 0& -0.01& 0.02\\ 0& 0& -0.02\end{array}\right)$(7.10)

where the diagonal element says how much the mass changes in that grading range and the off-diagonal terms say how much of a different size range is added to it. For instance, for the medium (i.e. i = 2$i=2$) size range

where $A_{2,2}=-0.01$ indicates that 1% of the mass in the medium size range is removed from it and $A_{3,2}=0.02$ indicates that 2% of the mass in the large size range is added to the medium size range.

Note that for mass conservation each column of the marginal transition matrix must add to zero. If mass is lost (e.g. by dissolution of particles), then the column(s) will sum to less than zero. Also note that none of the diagonal terms can be less than −1 because otherwise the equation is transforming more mass in the layer than actually exists in that layer. In this latter case a smaller timestep is required so that all the elements of the marginal transition matrix are smaller.

Finally it should be noted that, within the physical constraints above, there is considerable flexibility for the contents of the matrix, and therefore what processes can be simulated. The entries below the diagonal will normally all be zero because otherwise the process being represented by the matrix will be making larger particles from smaller particles (albeit if you are modelling soil aggregation or particle cementation this may be entirely reasonable). In the example in Equation (7.10) particles changed only to particles in the next size class down, which may not be the case, for instance, for particles that fragment into many smaller particles, or for spalling where there is a single large particle and many smaller particles resulting from the weathering. We will discuss these cases later.

While not presented this way the work of Salvador-Blanes et al. (2007) can be cast into this matrix formulation. Salvador-Blanes modified the approach of Legros and Pedro (Reference Legros and Pedro1985) to include the modelling of the breakdown of particles less than 2 mm in diameter and used Legros’s method whereby a particle was transformed into a particle in the next size class down (he defined that as 2 μm smaller). Both Legros and Salvador-Blanes used 1,000 size fractions for particles less than 2 mm, and this is easily transformed into Equation (7.10) (the A matrix will be 1, 000 × 1, 000$1,000\times 1,000$) with only the diagonal terms and that entry directly above the diagonal term being nonzero, and the values being the rates at which one size fraction is transformed to the next smaller size fraction per timestep.

The recasting of Salvador-Blanes’s approach highlights an important but implied aspect of the matrix methodology. The matrix A is not describing how an individual particle is breaking down, but the aggregate result of the breakdown of the many particles in a size fraction. It is not possible to move a single particle to the next size class smaller and maintain mass conservation without generating some other smaller particles. However, the matrix method describes the aggregate of a given mass of particles, and it is possible to have many particles transform to the next size class lower without generating an array of fine particles.

7.3.2 Interactions between Soil Layers

We now describe how one layer within the soil profile interacts with another layer within the soil profile. Let us consider two layers j and k. We can construct a matrix equation that describes how layer j changes layer k in one timestep:

${\underset{¯}{g}}_{k,t+1}=L_{k,j}{\underset{¯}{g}}_{j,t}$(7.12)

This equation says how much mass in layer j in each of the grading size fractions is added to or subtracted from layer k in one timestep. Equation (7.12) models only how material is moved between layers and not weathering. We have used the matrix notation L here to distinguish this movement between layers from the weathering process matrices (that transform the particle size distribution within a single layer). For mass conservation in layer k, if mass is added to layer k from layer j, then an equal amount of mass must be removed from layer k and moved to another layer. The combination of the weathering and movement can be expressed in matrix form if we construct a grading vector that merges all the grading vectors for each of the individual layers, and the marginal transition matrix is constructed from the layer transition matrices and interlayer movement transition matrices so that

${\left(\begin{array}{c}{\underset{¯}{g}}_1\\ {\underset{¯}{g}}_2\\ ⋮\\ {\underset{¯}{g}}_n\\ {\underset{¯}{g}}_{\infty }\end{array}\right)}_{t+1}={\left(\begin{array}{c}{\underset{¯}{g}}_1\\ {\underset{¯}{g}}_2\\ ⋮\\ {\underset{¯}{g}}_n\\ {\underset{¯}{g}}_{\infty }\end{array}\right)}_t+\left(\begin{array}{ccccc}{A}_{1,1}& L_{1,2}& \cdots & L_{1,n}& L_{1,\infty }\\ L_{2,1}& A_{2,2}& \cdots & L_{2,n}& L_{2,\infty }\\ ⋮& ⋮& \ddots & ⋮& ⋮\\ L_{n,1}& L_{n,2}& \cdots & A_{n,n}& L_{n,\infty }\\ \left(0\right)& \left(0\right)& \cdots & \left(0\right)& A_{\infty ,\infty }\end{array}\right){\left(\begin{array}{c}{\underset{¯}{g}}_1\\ {\underset{¯}{g}}_2\\ ⋮\\ {\underset{¯}{g}}_n\\ {\underset{¯}{g}}_{\infty }\end{array}\right)}_t$(7.13)

where the state vector is the gradings for all the layers from the surface armouring layer (subscript 1 indicates the surface armouring layer), through the profile layers and including the semi-infinite underlying layer (subscript ∞$\infty$). The state vector is thus a vector of length m(n + 1)$m\left(n+1\right)$ and the matrix is of dimension (m(n + 1)) × (m(n + 1)).$\left(m\left(n+1\right)\right)\times \left(m\left(n+1\right)\right).$ The notation [0]$\left(0\right)$ indicates a matrix that is m × m$m\times m$ and where all matrix entries are zeros. Equation (7.13) shows that all layers can interact with all other layers both above and below, including the semi-infinite underlying layer. Looking at the bottom row of the matrix, the semi-infinite layer can change through time (as a result of the matrix $A_{\infty ,\infty }$), but it cannot be influenced by any of the overlying layers (since all the entries are [0]$\left(0\right)$). Equation (7.13) can be written in a more compact form

${\underset{═}{g}}_{t+1}=\left(I+B\right){\underset{═}{g}}_t$(7.14)

where the double underbar notation distinguishes the grading vector for an individual layer (one underbar) from the vector for the gradings for all the layers (two underbars). The latter is sometimes called a supervector in the modelling literature (because it is a vector of vectors).

The construction of matrix B appears at face value to be rather daunting simply because of its size. However, Equation (7.13) is the most general statement of the problem, and in many cases simplifications are possible. For example,

• If interactions between layers occur only between adjacent layers: In this case for layer j, the only matrices that are nonzero are the matrices $A_{j,j}$ for the weathering within that layer, $L_{j,j-1}$ and $L_{j,j+1}$ which describe how the layer above and below, respectively change layer j, and $L_{j-1,j}$ and $L_{j+1,j}$ which describe how layer j changes the layer above and below, respectively:

  $\left(\begin{array}{ccccccc}{A}_{1,1}& L_{1,2}& \left(0\right)& \cdots & \left(0\right)& \left(0\right)& \left(0\right)\\ L_{2,1}& A_{2,2}& L_{2,3}& \cdots & \left(0\right)& \left(0\right)& \left(0\right)\\ \left(0\right)& L_{3,2}& A_{3,3}& \cdots & \left(0\right)& \left(0\right)& \left(0\right)\\ ⋮& ⋮& ⋮& \ddots & ⋮& ⋮& ⋮\\ \left(0\right)& \left(0\right)& \left(0\right)& \cdots & A_{n-1,n-1}& L_{n-1,n}& \left(0\right)\\ \left(0\right)& \left(0\right)& \left(0\right)& \cdots & L_{n,n-1}& A_{n,n}& L_{n,\infty }\\ \left(0\right)& \left(0\right)& \left(0\right)& \cdots & \left(0\right)& \left(0\right)& A_{\infty ,\infty }\end{array}\right)$(7.15)

• If in addition to interactions only between adjacent layers, there is no change in the grading over time within a layer: This might happen when there is only mixing of the soil (e.g. bioturbation) and no breakdown of the mineral matter. The matrix in Equation (7.15) simplifies even further to

  $\left(\begin{array}{ccccccc}\left(0\right)& L_{1,2}& \left(0\right)& \cdots & \left(0\right)& \left(0\right)& \left(0\right)\\ L_{2,1}& \left(0\right)& L_{2,3}& \cdots & \left(0\right)& \left(0\right)& \left(0\right)\\ \left(0\right)& L_{3,2}& \left(0\right)& \cdots & \left(0\right)& \left(0\right)& \left(0\right)\\ ⋮& ⋮& ⋮& \ddots & ⋮& ⋮& ⋮\\ \left(0\right)& \left(0\right)& \left(0\right)& \cdots & \left(0\right)& L_{n-1,n}& \left(0\right)\\ \left(0\right)& \left(0\right)& \left(0\right)& \cdots & L_{n,n-1}& \left(0\right)& L_{n,\infty }\\ \left(0\right)& \left(0\right)& \left(0\right)& \cdots & \left(0\right)& \left(0\right)& \left(0\right)\end{array}\right)$(7.16)

• Adjusting the layers in response to erosion at the surface: When material is removed from the surface layer, then material to balance that lost to erosion must be removed from the layer below to make sure the mass in the layer does not change. This cascades down through all layers in the profile:

  $\left(\begin{array}{ccccccc}-\frac{E}{d_1}{A}_{1,1}& \frac{E}{d_1}I& \left(0\right)& \cdots & \left(0\right)& \left(0\right)& \left(0\right)\\ \left(0\right)& -\frac{E}{d_2}I& \frac{E}{d_2}I& \cdots & \left(0\right)& \left(0\right)& \left(0\right)\\ \left(0\right)& \left(0\right)& -\frac{E}{d_3}I& \cdots & \left(0\right)& \left(0\right)& \left(0\right)\\ ⋮& ⋮& ⋮& \ddots & ⋮& ⋮& ⋮\\ \left(0\right)& \left(0\right)& \left(0\right)& \cdots & -\frac{E}{d_{n-1}}I& \frac{E}{d_{n-1}}I& \left(0\right)\\ \left(0\right)& \left(0\right)& \left(0\right)& \cdots & \left(0\right)& -\frac{E}{d_n}I& \frac{E}{d_n}I\\ \left(0\right)& \left(0\right)& \left(0\right)& \cdots & \left(0\right)& \left(0\right)& \left(0\right)\end{array}\right)$(7.17)

where E is the erosion in one timestep in depth units, d_i is the depth of layer i (which converts sediment mass due to erosion E into a proportion of the layer mass) and I is a k × k identity matrix. The matrix A is the armouring transition matrix for the surface layer and determines the size selectivity of the sediment entrainment due to erosion. As an aside, this is the first time that bulk density of the soil appears in the matrix methodology. Bulk density is the conversion factor between depth of soil eroded and mass of soil eroded. If soil erosion is expressed in mass rather than depth units, then the bulk density conversion is not required.

7.3.3 Constructing the A Matrix for Weathering

We will now discuss how the matrix A is formulated to model weathering. Unless otherwise stated, the discussion in this section is for one layer, and we drop the layer subscript j for clarity of discussion. We will first discuss mass conservative weathering and then generalise the discussion to non-mass conservative weathering afterwards.

The basis of the conceptualisation of weathering is that particles are spread uniformly within each size range within which they fall (Figure 7.2). Thus some particles will be at the lower boundary of the size, some in the middle and some at the top of the size range. We conceptualise the fragmentation process as a parent particle breaking into a number of daughter particles (Figure 7.3). Figure 7.3 shows that depending on the fragmentation mechanism there may be a range of different types of daughter particles created. If we assume that the density distribution within a single particle size fraction is constant, then mass conservation implies volume conservation. In this case the total volume of the daughter particles is the same as the volume of the parent particle.

Figure 7.2: How soil grading is conceptualised in mARM and SSSPAM as uniformly distributed within size classes.

Data are the Ranger Mine grading used in Willgoose and Sharmeen (Reference Willgoose and Sharmeen2006) and Cohen et al. (Reference Cohen, Willgoose and Hancock2009, Reference Cohen, Willgoose and Hancock2010).

(from Sharmeen and Willgoose, Reference Sharmeen and Willgoose2006)

Figure 7.3: Schematic of fracturing models for rock particles (a) breaking of small particles from a rind and (b) body fragmentation.

We now extend our discussion from what happens to a single particle to what happens to all the particles within a size class range. Figure 7.4 shows what will happen to the larger size range when the particles break into two particles of different diameter, if they all follow the rules for breaking of a single particle (i.e. all particles break with exactly the same geometry). In the discussions that follow we assume that all particles are spherical. Initially we will consider the case where a parent particle breaks into two daughter particles. Figure 7.4 shows how the geometry of a single particle breaking allows us to map the parent particle size grade to the smaller dimensions of the two daughter particles’ size grades. Likewise we can map the largest size from the parent size gradings to the largest size of the daughter size grading. From the geometry of breaking of the single particle we know what the proportions of the mass of daughter particles should be in each of the two daughter particle gradings ranges. Figure 7.4 shows that if we take one parent size fraction, then the daughter particle’s size distribution will not typically neatly fit into a single size fraction, but will need to be interpolated across a number of size fractions. These proportions will be a function of the diameters of the daughter particles and of the size grading fractions adopted by the user in the modelling. Finally, if the distribution between the lower and upper size of the parent particles is uniform, then the distribution of the daughter particles is also uniform between the lower and upper size grading.

Figure 7.4: Conceptualisation of how size fractions in the grading are transformed in each time step of mARM.

(from Cohen et al., Reference Cohen, Willgoose and Hancock2009)

Using these assumptions it is relatively straightforward, though tedious, to construct a matrix to simulate weathering. The calculations in Figure 7.4 are done for each grading range and the results summed.

However, a conceptual simplification is possible. There is a combination of grading size fractions and weathering processes that leads to particularly simple results, and which allow us to derive some analytical results that provide useful insight into the time variation of grading under the action of weathering. We will use this to demonstrate how the weathering matrix works and derive some simple weathering results. We will then show that this simple model can be used to derive more complex models.

For this example we will define the fractions in the grading size distribution such that the lower diameter limit of each size fraction is $\left({0.5}^{\frac{1}{3}}\right)$ times the upper diameter limit for that same size fraction. If we take the maximum size limit as 2 mm, then this yields the limits for the size fractions as in Table 7.1. There is no special significance to the upper and lower limits of 2 mm and 0.125 mm in the table; they are simply used to show what this grading looks like. It is useful to examine this size grading because this grading fractionation is different from that normally used for soils analysis (e.g. phi grading) and thus will be unfamiliar to readers. The significance of the grading in this example is that if a spherical particle splits into two equal volume spherical daughter particles, then for conservation of volume all the particles in one size grading will fall exactly into the size grading of the next size fraction smaller. For example, for the grading class 12 in Table 7.1 a particle with the largest diameter of 2 mm will split into two particles of diameter 1.587 mm, and a particle with the smallest diameter of 1.587 mm will split into two particles of diameter 1.260 mm. These upper and lower size limits define the size range for the next smaller Grading Class 11. This is true for all Grading Class ranges in Table 7.1.

Accordingly, to assemble the weathering matrix with this grading when a parent particle fractures into two equally sized daughter particles is straightforward because the daughter particles only ever fall into the size fraction below, and this size fraction below can receive only daughter particles that fracture from the size fraction directly above. Thus if we say that the proportion of the particles that fracture for one timestep is α,$\alpha ,$ then the weathering matrix is

$A=\left(\begin{array}{ccccccc}0& \alpha & 0& \cdots & 0& 0& 0\\ 0& -\alpha & \alpha & \cdots & 0& 0& 0\\ 0& 0& -\alpha & \cdots & 0& 0& 0\\ ⋮& ⋮& ⋮& \ddots & ⋮& ⋮& ⋮\\ 0& 0& 0& \cdots & -\alpha & \alpha & 0\\ 0& 0& 0& \cdots & 0& -\alpha & \alpha \\ 0& 0& 0& \cdots & 0& 0& -\alpha \end{array}\right)$(7.18)

The diagonal and off-diagonal α terms mean that α of the mass in the larger size grading is all added to the next smaller size grading. The diagonal element for the smallest fraction is zero because we assume that particles this size cannot weather any smaller. Equation (7.18) is mass conservative, so all the columns of the matrix add to 0.

It is possible to calculate the time variation of the mean diameter of the particle size distribution purely from the I + A$I+A$ matrix using Equation (7.18). With the exception of the smallest size grading, the mass g_i$g_i$ that is in any given size fraction i will change after one timestep so that (1 − α)$\left(1-\alpha \right)$ remains in size fraction i while α$\alpha$ will now be in the next size fraction smaller i − 1,$i-1,$ which in the case of Equation (7.18) is $\left({0.5}^{\frac{1}{3}}\right)$ smaller than the size fraction above, and from the way the size fractions are defined in Table 7.1 that is true for all size fractions i. For the smallest size fraction there is no change. If we consider the case where all the mass is initially concentrated in the largest size fractions and there is an insignificant mass in the smallest fraction (this can be achieved by defining the smallest size fraction to be very small), then in one timestep (time changing from t to t + 1$t+1$) the mean diameter of the soil changes from

$d_{\text{mean},t+1}=\left(\left(1-\alpha \right)+\alpha {0.5}^{\frac{1}{3}}\right)d_{\text{mean},t}=\left(1-0.2063\alpha \right)d_{\text{mean},t}$(7.19)

which allows us to derive an equation for the evolution of the mean diameter of the soil T timesteps into the future:

which is a semi-log linear relationship between diameter and time, where the slope on a semi-log plot is −0.2063α$-0.2063\alpha$. This can also be expressed as an exponential so that

Figure 7.5 shows this curve for a number of values of α$\alpha$ and with a starting mean diameter of 3.5 mm.

Figure 7.5: Numerical solution of the time evolution of the mean of the particle size distribution for a range of weathering rates. The parameter ‘rate’ is rate = 0.2063α$\text{rate}=0.2063\alpha$ in Equation (7.19). The levelling off on the right-hand side of each graph is because at this stage all soil is in the finest size fraction in the modelled soil grading.

The exact details of Equations (7.19) to (7.21) are a function of the specific assumptions in Equations (7.18), but it is straightforward to extend this analysis to any distribution of daughter products, and the only part of Equation (7.19) that will change is the number inside the parentheses; the relationship itself will still be log-log linear. This is true provided only that (1) the definition of the size fractions is defined as in Table 7.1 for two daughter particles, (2) the fracture model is independent of the diameter of the particle so that the daughter particles are always the same size relative to the parent particle and that this fracture model is independent of the diameter of the parent particle and (3) the smallest size fraction is small enough that only a minor amount of the mass occurs within the smallest fraction (otherwise the curves will start to level out as the particles age as seen in Figure 7.5 on the right-hand side).

The assumption that particles split into two equally sized particles leading to Equation (7.21) is rather restrictive, but this result can be generalised as follows. Consider first the case where N equally sized particles are generated at each time step instead of the two in the example above. This means that the diameter of the daughter particles is now $N^{\frac{-1}{3}}$ the diameter of the parent particle. In the same way that we developed a grading range definition for two daughter particles fracturing in Table 7.1, we can do this in general for all values of N. In Table 7.1 fracturing for three equal daughters is shown. Note that the total size range covered by the 12 size fractions in Table 7.1 is different, with the three daughters column going down to 0.025 mm, while the two daughters range goes down to only 0.125 mm. Thus if you wish to cover a specific total size range (e.g. for a given soil), the two daughters size range requires more size fractions. If the proportion of a given size fraction that breaks down each timestep is α$\alpha$ as before, then Equation (7.18) is still applicable except now the definition of the size grading ranges is different and the matrix will be of a different dimension. The equivalent three daughters result to Equation (7.19) is

$d_{\text{mean},t+1}=\left(\left(1-\alpha \right)+{0.333}^{\frac{1}{3}}\alpha \right)d_{\text{mean},t}=\left(1-0.3066\alpha \right)d_{\text{mean},t}$(7.22)

and Equation (7.20) becomes

In general for N daughter particles

$d_{\text{mean},t=T}={\left(1-\left(1-N^{\frac{-1}{3}}\right)\alpha \right)}^T{d}_{\text{mean},t=0}$(7.24)

so it is clear that if particles break into a larger number of smaller particles at each timestep, then the soil will become finer at a faster rate, and Equation (7.24) defines what that speedup will be.

In the discussion that follows we further generalise the fracture geometries considered. To simplify the discussion of geometry we introduce a shorthand for fragmentation geometry, which indicates how, on average, the volume of particles will be distributed after a fragmentation event assuming the size grading definitions in Table 7.1. The general form is pX-AAA-BBB-CCC-DDD where X is the number of daughter particles the grading fractions have been defined for (i.e. X = 2$X=2$ for the two particle grading in Table 7.1, X = 3$X=3$ for the three particle grading), AAA is the percentage of volume that remains in the parent size grading fraction after fragmentation, BBB is the proportion in the next size grading smaller, CCC in the next size down again and so on. In this notation where all of the particles split into two with two equally sized particles, then the shorthand is p2-0–100. Some geometry examples are listed in Table 7.2. Some simple examples of fracture geometries and their fragmentation notation are also shown. While all the examples in Table 7.1 have no mass left in the parent grading after fragmentation (they all have a leading zero), at the end of this section an important, but more complex, fragmentation geometry that leads to a nonzero percentage in the parent grading will be discussed.

Instead of two equally sized daughter particles, let us now break the particle into three daughters where one daughter is half the volume of the original and the other two daughters are a quarter of the volume of the original (i.e. p2-0–50-50). Mass conservation still applies, but what now happens is that half the weathered volume goes into the next class down from the parent, while half the weathered volume goes into the class 2 classes below the parent so the A matrix is

$A=\left(\begin{array}{ccccccc}0& \alpha & \frac{\alpha }{2}& \cdots & 0& 0& 0\\ 0& -\alpha & \frac{\alpha }{2}& \cdots & 0& 0& 0\\ 0& 0& -\alpha & \cdots & 0& 0& 0\\ ⋮& ⋮& ⋮& \ddots & ⋮& ⋮& ⋮\\ 0& 0& 0& \cdots & -\alpha & \frac{\alpha }{2}& \frac{\alpha }{2}\\ 0& 0& 0& \cdots & 0& -\alpha & \frac{\alpha }{2}\\ 0& 0& 0& \cdots & 0& 0& -\alpha \end{array}\right)$(7.25)

and Equation (7.19) becomes

$\begin{array}{cc}{d}_{\text{mean},t+1}& =\left(\left(1-\alpha \right)+\frac{\alpha }{2}{0.5}^{\frac{1}{3}}+\frac{\alpha }{2}{0.25}^{\frac{1}{3}}\right)d_{\text{mean},t}\\ & =\left(1-0.2882\alpha \right)d_{\text{mean},t}\end{array}$(7.26)

As expected this case weathers faster than the two equal daughters case. In a similar fashion you could consider that this three daughters case can be extended where one of the two smaller particles itself breaks into two, so the daughter particles are 1 × (1/2)$1\times \left(1/2\right)$ volume, 1 × (1/4)$1\times \left(1/4\right)$ volume and 2 × (1/8)$2\times \left(1/8\right)$ volume and so on.

Similarly the three equal daughters particle fracturing could be extended by having one of the daughters break into three equal volume smaller particles, so that the daughters are two particles with 1/3 the parent volume (i.e. 2/3 of the parent volume) and three particles with 1/9 the parent volume (i.e. 1/3 of the parent volume), yielding p3-67.7-33.3. In this case Equation (7.25) becomes

$A=\left(\begin{array}{ccccccc}0& \alpha & \frac{\alpha }{3}& \cdots & 0& 0& 0\\ 0& -\alpha & \frac{2\alpha }{3}& \cdots & 0& 0& 0\\ 0& 0& -\alpha & \cdots & 0& 0& 0\\ ⋮& ⋮& ⋮& \ddots & ⋮& ⋮& ⋮\\ 0& 0& 0& \cdots & -\alpha & \frac{2\alpha }{3}& \frac{\alpha }{3}\\ 0& 0& 0& \cdots & 0& -\alpha & \frac{2\alpha }{3}\\ 0& 0& 0& \cdots & 0& 0& -\alpha \end{array}\right)$(7.27)

and Equation (7.26) becomes

$\begin{array}{cc}{d}_{\text{mean},t+1}& =\left(\left(1-\alpha \right)+\frac{2\alpha }{3}{0.333}^{\frac{1}{3}}+\frac{\alpha }{3}{0.333}^{\frac{2}{3}}\right)d_{\text{mean},t}\\ & =\left(1-0.3775\alpha \right)d_{\text{mean},t}\end{array}$(7.28)

If two of the three daughters break into three (i.e. p3-33.3–66.7), then

$\begin{array}{cc}{d}_{\text{mean},t+1}& =\left(\left(1-\alpha \right)+\frac{\alpha }{3}{0.333}^{\frac{1}{3}}+\frac{2\alpha }{3}{0.333}^{\frac{2}{3}}\right)d_{\text{mean},t}\\ & =\left(1-0.4484\alpha \right)d_{\text{mean},t}\end{array}$(7.29)

and as for the two daughters cases, both of these cases weather faster than the case where all the daughters are of equal sizes, and when more fine particles are generated (i.e. Equation (7.29) versus (7.28), or Equation (7.28) versus (7.26)) then the soil weathers faster even when the rate of breakdown per timestep is unchanged.

It should be clear from these examples that the combinations of sizes of daughter particles possible are quite extensive. The only constraint is a geometry constraint driven by the grading size definition adopted from Table 7.1. If the two daughter grading range in Table 7.1 is adopted, then particles must break into particles where the volumes of the daughters are related to the parent particle volume by integer powers of 1/2. Likewise for the three daughter grading range in Table 7.1, then all the daughters must have volumes that are integer powers of 1/3. This integer power constraint is so that all the daughters fall entirely and only into one of the grading size ranges and do not span two grading size ranges.

This section has shown how to construct the weathering A matrix using a conceptualisation of the fracturing of individual parent particles into a range of daughter particles. Thus everything in this section has been strongly physically based on fracturing mechanisms that are, in principle, observable in the laboratory or the field. This presentation was intentional because it highlights the link between fracturing of individual particles and the rather less tangible mathematics of the A matrix.

We will deviate from this philosophy for a moment to generalise Equation (7.18), which will then allow us to provide a general analytical solution to the change in grading over time. When constructing the A matrix, to ensure mass conservation we need to ensure only that (1) the summation of each column should be zero, (2) the diagonal term is between −1 and 0 and (3) the off-diagonal terms are greater than or equal to 0. Thus if we can construct a generic column for matrix A, that is the same for all grading fractions

$\text{column}i\text{in}A=\left(\begin{array}{c}{a}_{i,1}\\ ⋮\\ a_{i,i-2}\\ a_{i,i-1}\\ a_{i,i}\\ 0\\ 0\\ ⋮\\ 0\end{array}\right)$(7.30)

we can then write Equation (7.19) as (using the two daughter size grading)

$d_{\text{mean},t+1}=\left(\left(1-a_{i,i}\right)+a_{i,i-1}{0.5}^{\frac{1}{3}}+a_{i,i-2}{0.5}^{\frac{2}{3}}+\dots +a_{i,1}{0.5}^{\frac{\left(i-1\right)}{3}}\right)d_{\text{mean},t}$(7.31)

and if the values of a_{i, j}$a_{i,j}$ are small when j is small (i.e. large parent particles do not generate many very small daughter particles), then we can write

$d_{\text{mean},t+1}\cong \overline{a}{d}_{\text{mean},t}$(7.32)

where $\overline{a}$ is a constant reflecting the structure of A so that we can write a general equation for how the mean diameter of the soil will change over time

$d_{\text{mean},t=T}={\left(\overline{a}\right)}^T{d}_{\text{mean},t=0}$(7.33)

which shows that the mean of the soil grading is an exponential function with time of the initial grading irrespective of the structure of A. Note that the restrictions on the form of A are relatively modest (just that weathering cannot generate a large mass proportion of small particles), so this is quite a powerful and general conclusion. Note also that the use of the two daughter size fraction definition in Table 7.1 is a convenience rather than a necessity, as we will see in the next section. The fragmentation notation for this case is p2-(1-a_i,i)-a_i,i-1- ^… -a_i,1.

The discussion above presumes that all particles fracture with exactly the same geometry. This is not essential. One simple extension is to allow particles to break with two possible geometries at different rates. The A matrix for the combined process is then

$A=w_1{A}_1+w_2{A}_2$(7.34)

where the subscripts are the weathering rate w and weathering matrix A for each of the two processes. This can be generalised to as many processes as required.

Figure 7.6 extends the discussion by examining a number of fracture geometries, and shows how (1) the different fracture geometries yield different rates of evolution of the mean grading and (2) all the geometries evolve as a semi-log of time as in Equation (7.33). The meaning of the coding used for describing fracture geometry is explained in the figure caption, but is based on the fracture geometry where particles are 1/2, 1/4, 1/8 and so on of the volume of the parent particle.

Figure 7.6: (a) The average particle size distribution of the daughter particles resulting from random fragmentation (P2-stoch), d_parent is the diameter of the parent particle, (b) Numerical solution of the time evolution of the mean of the particle size distribution for the same weathering rate (as used in Figure 7.5) and a range of fragmentation models. The weathering model p2-0–50-50 is the same as ‘rate = 0.03’ in Figure 7.5.

One of the fracture geometries in Figure 7.6 is philosophically different from the ones described above. In all the fracture geometries above the geometry of the daughter particles is deterministically fixed as a proportion of the volume of the parent particle. The aggregated behaviour of all the particles is then just the sum of all the particles in that grading range. However, consider the case where a single particle breaks into two particles but where the fracture location in the particle is randomly distributed within the particle (e.g. a schist where the particle cleaves along the layering but where the cleavage plane is randomly located within the particle). The probability distribution of the volume of daughter particles derived from particles in a specified particle range is shown in Figure 7.6a. It is convenient here to define the x axis as the volume of particles rather than their diameter because the probability distribution of the daughter particles is then particularly simple. If a large number of particles fracture, then the probability distribution function gives the mean of the volume distribution of the daughter particles, from which A can be easily derived. In fragmentation notation it is p2-30.7-34.6-17.3-17.4, if we lump the volume of all particles more than three grading fractions smaller into the third grading fraction (i.e. the 17.4). Other than the proportion remaining in the parent grading range, the relativities between the daughter size fractions are the same as the scaling model p2-0–50-25-12.5-12.6.

Figure 7.6 shows that those fractions that generate a greater proportion of volume in smaller size fractions fine faster, which is consistent with Equation (7.24). Figure 7.7 shows the particle size distribution resulting from using a number of different fragmentation models. Initially all particles were 4.0 mm in diameter, and Figure 7.7 shows the particle size distribution when the fragmentation models have weathered the particles to a median diameter 0.1 mm. Two conclusions of this figure are that (1) different fragmentation models yield slightly different particle size distributions even though the median diameter of the particle size distribution is the same and (2) while the different fragmentation models may yield finer particle size distributions at different rates, they yield particle size distributions that are only slightly different (e.g. the p2-0–100 and the p2-stoch random fracture models are indistinguishable in Figure 7.7, but p2-0–100 fines a factor of two times faster in Figure 7.6).

Figure 7.7: These are the same simulations as in Figure 7.6 (i.e. same weathering rate, different fragmentation models) but showing the particle size distribution when the d₅₀ of the all fragmentation models is equal to 0.1 mm. Note that Figure 7.6 shows that the age of each of the distributions will be different since the different fragmentation models generate a different rate of change in the d₅₀ for the same weathering rate.

7.3.4 Transforming the Size Grading Definition

The definition of the size fractions necessary to derive Equations (7.18) to (7.29) is rather restrictive, and is not necessarily consistent with size gradings used in soil science practice (e.g. sieve sizes). Accordingly, we will now show how to generalise the definition of diameters defining the size fractions while still maintaining the simplicity of the equation and the semi-log relationship for mean size with time. It is possible to transform the mass fractions from one size fraction definition to another size fraction definition. Consider $\underset{¯}{g}$ to be the grading by one grading class size definition and $\underset{¯}{h}$ to be the desired mass fraction equivalent to $\underset{¯}{g}$ in another grading class size definition. For generality let us consider that there are n size classes in $\underset{¯}{g}$ and m size classes in $\underset{¯}{h}$. We can transform $\underset{¯}{g}$ into $\underset{¯}{h}$ with the matrix equation

$\underset{¯}{h}=T\underset{¯}{g}$(7.35)

where T is a n × m matrix that distributes the mass in the n classes in $\underset{¯}{g}$ into the m size classes in $\underset{¯}{h}$.

We can now transform the weathering (or any other process for that matter) matrix equation that is originally

${\underset{¯}{g}}_{t+1}=\left(I+A\right){\underset{¯}{g}}_t$(7.36)

By substituting in the inverse of Equation (7.35) we get the transition matrix in the new grading

$\begin{array}{c}{T}^{-1}{\underset{¯}{h}}_{t+1}=\left(I+A\right)T^{-1}{\underset{¯}{h}}_t\\ {\underset{¯}{h}}_{t+1}=T\left(I+A\right)T^{-1}{\underset{¯}{h}}_t=\left(I+\hat{A}\right){\underset{¯}{h}}_t\end{array}$(7.37)

so that in the new grading definition

$I+\hat{A}=T\left(I+A\right)T^{-1}=I+{\mathrm{TAT}}^{-1}$(7.38)

noting that ${\mathrm{TT}}^{-1}=I.$

This shows that we can populate the A matrix using any size class definition we like and then later transform it into any other size class definition that is convenient. One major consequence of this is that the semi-log evolution of the mean diameter we have derived in the previous section is not unique to the size class distribution used to derive it, and that even if the class definition is different (provided the physics is the same), then the semi-log relationship with time will still be observed.

7.3.5 Constructing the A Matrix for Pedoturbation

Pedoturbation is any process that disturbs the horizonation of soils. Some pedoturbation processes assist in the differentiation of soils into layers, while others mix the soils, tending to destroy the layers (Hupy and Schaetzl, Reference Hupy and Schaetzl2006). Bioturbation is a pedoturbation process resulting from biological activity and is a significant process in soil profile development (Hole, Reference Hole1981). One of the main impacts of bioturbation is to physically mix soil from top to bottom (Heimsath et al., 2002; Ahr et al., Reference Ahr, Nordt and Forman2013), and transport material from the surface into the profile (Astete et al., Reference Astete, Constant, Thibodeaux, Seals and Delim2015). Soil fauna such as ants and termites bring fine material from deep in the soil to the surface (Lobry de Bruyn and Conacher, Reference Lobry de Bruyn and Conacher1990), earthworms likewise mix the profile (Van Hooff, Reference van Hooff1983; McKenzie et al., Reference McKenzie and Dexter1993), while tunnelling animals such as gophers (Huntley and Inouye, Reference Huntly and Inouye1988; Gabet, Reference Gabet2000) and wombats (Paton et al., Reference Paton, Humphreys and Mitchell1995; Field and Anderson, Reference Field and Anderson2003) may also transfer material laterally (via their horizontal tunnels) as well as from within the profile to the surface. Flora such as trees mix the soil through tree throw (Gabet et al., Reference Gabet and Dunne2003; Field and Anderson, Reference Field and Anderson2003). In this section we will consider only vertical transport of material. Horizontal movement will be discussed in Chapter 9.

Many of these processes leave macropores within the profile when material is removed from that layer, which means that the porosity is increased and the bulk density of the soil layers is reduced. If the mass in each soil layer in the model is the same and bulk density decreases, then the layers will dilate (or expand), so that the soil surface rises. For instance, Wilkinson et al. (Reference Wilkinson, Richards and Humphreys2009) provides data for bioturbation rates that decline with depth roughly exponentially from the surface, which they attribute primarily to ants and worms. They also measured bulk density and found increases with depth (1.1 g/cm³ near to the surface to 1.8 g/cm³ at a depth of 0.8 m) mirroring the decline in the rate of bioturbation and number of macropores generated. Shull (Reference Shull2001) presents a very similar matrix model for bioturbation in sea bottom sediments.

One unusual form of pedoturbation was the soil mixing resulting from artillery shells in France during World War I (‘bombturbation’, Hupy and Schaetzl, Reference Hupy and Schaetzl2006, Reference Hupy and Schaetzl2008), and these and other warfare sites provide unique field experiments in pedogenesis (Certini et al., Reference Certini, Scalenghe and Woods2013).

7.3.5.1 Simple Biological Mixing

A simple model for bioturbation mixes the soil vertically randomly, and this can be modelled as a vertical diffusive process. The idea is that microfauna (e.g. earthworms) mix the soil vertically. This process could be modelled by a Fickian diffusion-type equation. A simple method, which asymptotes to Fickian diffusion (i.e. $Q_s=D\frac{\text{∂}C}{\text{∂}x}$ where Q_s$Q_s$ is the rate of vertical transport per unit area and time, $\frac{\text{∂}C}{\text{∂}x}$ is the gradient of the property C being mixed and D is the diffusivity), is to exchange material with the layers above and below (Williams, Reference Williams2006):

$B=\left(\begin{array}{cccccc}-R_{1}I& R_{2}I& \left(0\right)& \left(0\right)& \left(0\right)& \cdots \\ R_{1}I& -2R_{2}I& R_{3}I& \left(0\right)& \left(0\right)& \cdots \\ \left(0\right)& R_{2}I& -2R_{3}I& R_{4}I& \left(0\right)& \cdots \\ \left(0\right)& \left(0\right)& R_{3}I& -2R_{4}I& \ddots & \cdots \\ \left(0\right)& \left(0\right)& \left(0\right)& R_{4}I& \ddots & \ddots \\ ⋮& ⋮& ⋮& ⋮& \ddots & \ddots \end{array}\right)$(7.39)

where R_i is the rate of bioturbation mixing at the depth of layer i. Equation (7.39) assumes that all layers have the same thickness. Following Wilkinson et al. (Reference Wilkinson, Richards and Humphreys2009) the bioturbation rate R_i declines with depth approximately exponentially. The identity matrix I means that the mixing is not size selective. For constant thickness T layers, this converges to Fickian diffusion asymptotically with time

$R_i=D\frac{\Delta t}{T^2}$(7.40)

If the layers are of different thickness, then the ith column of tridiagonal elements changes to

$\left(\begin{array}{c}\frac{2D\Delta t}{{\left(\Delta z_1\right)}^2\left(1+\frac{1}{\alpha }\right)}\\ -\frac{2D\left(1+\alpha \right)\Delta t}{{\left(\Delta z_1\right)}^2\left(1+\frac{1}{\alpha }\right)}\\ \frac{2\mathrm{D\alpha }\Delta t}{{\left(\Delta z_1\right)}^2\left(1+\frac{1}{\alpha }\right)}\end{array}\right)$(7.41)

where (see definitions in Figure 7.8) Δz₁$\Delta z_1$ is the distance between the nodes i and i − 1,$i-1,$ Δz₂$\Delta z_2$ is the distance between nodes i + 1$i+1$ and i, and $\alpha =\frac{\Delta z_1}{\Delta z_2}$. The relationship in Equation (7.41) is derived from a central finite difference approximation to Fickian diffusion. When the thicknesses of the three layers are the same (i.e. Δz₁ = T_i$\Delta z_1=T_i$ and α = 1$\alpha =1$) Equation (7.41) reduces to Equation (7.40).

Figure 7.8: Definitions of layer thickness and node discretisation for Equation (7.41).

7.3.5.2 Tree Throw

Tree throw brings to the surface all the soil in the root ball of the tree leaving behind a large hole in the ground. The soil falls off the roots over time, and the mixed soil in the root ball refills the hole in the ground. In addition to the mixing, there is potentially an impact on the porosity of the soil with the soil becoming more disordered. Tree throw is an event: a tree falls, and the consequences happen immediately. It is not a slow ongoing process, though the cumulative impact of many tree throw events is an ongoing process. Tree throw can be modelled as an event if it is assumed that the refilling of the hole happens quickly relative to the timestep of the modelling. At a given time a tree falls at a particular location or it doesn’t. Thus we do not use Equation (7.14) as the model of what happens from t to t + 1,$t+1,$ but we can use it to model what happens before and after the tree fall event

${\underset{═}{g}}_{\text{after}}=\left(I+B\right){\underset{═}{g}}_{\text{before}}$(7.42)

We simply need to know, when a tree fall event occurs, what the B$\text{B}$ matrix looks like for a single tree fall event, and apply Equation (7.42) for that time (and apply the equation repeatedly if there are many tree fall events at the same time). The mixing of the soil over the depth of the root zone is modelled using Equation (7.42) if we assume that the root ball rips up the top n layers so the top n layers are mixed together and the n + 1$n+1$ and lower layers are undisturbed. This is an averaging of the grading of the top n layers so that the B$\text{B}$ matrix can be expressed as

$B=\left(\begin{array}{cccccc}\left(\frac{m_1}{\sum _{k=1}^n{m}_k}-1\right)I& \cdots & \frac{m_1}{\sum _{k=1}^n{m}_k}I& \left(0\right)& \left(0\right)& \left(0\right)\\ ⋮& \ddots & ⋮& \left(0\right)& \left(0\right)& \left(0\right)\\ \frac{m_n}{\sum _{k=1}^n{m}_k}I& \cdots & \left(\frac{m_n}{\sum _{k=1}^n{m}_k}-1\right)I& \left(0\right)& \left(0\right)& \left(0\right)\\ \left(0\right)& \left(0\right)& \left(0\right)& \left(0\right)& \left(0\right)& \left(0\right)\\ \left(0\right)& \left(0\right)& \left(0\right)& \left(0\right)& \left(0\right)& \left(0\right)\\ \left(0\right)& \left(0\right)& \left(0\right)& \left(0\right)& \left(0\right)& \left(0\right)\end{array}\right)$(7.43)

where m_j is the total mass of layer j. The zero matrices for i greater and equal to n + 1$n+1$ indicate that the tree throw does not change the grading for layers n + 1$n+1$ and deeper. The term $\frac{m_j}{\sum _{k=1}^n{m}_k}$ is included to ensure that the mass in layer j is the same before and after the tree fall event. This is not critical but is convenient, because it preserves the thinner armour layer, which would otherwise become the same thickness as the underlying layers if the mass were distributed equally. The −1 in the diagonal term ensures that the grading in that layer is entirely replaced by the new grading to negate the effect of the identity matrix in Equation (7.42).

When modelling soil spatially it will often be the case that a single tree throw event disturbs only a proportion of the area associated with each computational node. If we call this proportion α,$\alpha ,$ then the equation for that node is a modification of Equation (7.42) so that

${\underset{═}{g}}_{\text{after}}=\left(\left(1-\alpha \right)I+\alpha B\right){\underset{═}{g}}_{\text{before}}$(7.44)

where the (1 − α)$\left(1-\alpha \right)$ term is the proportion of the nodal area that is unchanged by the tree throw event while the α$\alpha$ term is that proportion of the nodal area impacted by tree throw. The parameter α$\alpha$ might change from tree throw event to tree throw event reflecting the size of the tree (and thus the root ball) impacted. Likewise the depth of mixing might change between events to reflect tree throw for different size trees.

When modelling the temporal effect of tree throw (i.e. over a number of years rather than over a single tree throw event), there are two main alternative approaches.

The first approach is to model tree throw as occurring randomly in time and space (i.e. a Poisson random process) and modelling every tree throw event using the equations above. If tree throw has a probability of occurring of p per unit area per timestep, then the probability P of a tree throw event in a given timestep in a node of area A is P = pA.$P=\mathrm{pA}.$ The simplest way to simulate this is to cycle through every node at each timestep and randomly generate a uniformly distributed random number between 0 and 1 for each node. If the random number is greater than P, then a tree throw event occurs in that node, while if it is less than P, a tree throw event does not occur. If the area of a node is large, and/or the timestep is long, then it is possible for P to be greater than 1 (which indicates more than one tree throw event). In this case either one of the timestep or node area needs to be decreased.

The second approach calculates the average effect of tree throw over many events and models the mean behaviour rather than the effect of individual events. In this case the mean of the tree throw events in any given node is the rate of tree throw occurrence per unit time, which is P, so that the equation for the effect of tree throw per unit time, by modifying Equation (7.44), becomes

${\underset{═}{g}}_{t+1}=\left(\left(1-\mathrm{\alpha P}\right)I+\mathrm{\alpha P}B\right){\underset{═}{g}}_t$(7.45)

The advantage of the first method is that you can capture the stochastic effect (in both space and time) of tree throw, whereas the second method gives only the average effects. This stochastic effect may be important in the understanding of the spatial variability of the soil properties (e.g. Samonil et al., Reference Samonil, Kral and Hort2010; Gabet and Mudd, Reference Gabet and Mudd2010). The average of the first method gives the second method.

Because both these methods independently examine each node, the user can impose spatial variability on tree throw reflecting that tree throw is spatially clustered around the edges of forests and in openings within the forest, and that throw events are clustered in time during high-wind events. In this case the user can construct a wind event model, where wind is randomly generated across the domain, and based on the wind speeds impose a spatial pattern of tree throw occurrence (e.g. Espirito-Santo et al., Reference Espirito-Santo2014) and study its effect on the spatial variability of soil development (Finke et al., Reference Finke, Vanwalleghem, Opolot, Poesen and Deckers2013).

7.3.5.3 Termites and Ants

Termites and ants move fine soil material from the subsurface to the surface as part of nest construction (Paton et al., Reference Paton, Humphreys and Mitchell1995; Zaitlin and Hayashi, Reference Zaitlin and Hayashi2012). The surface expression of these nests is then, over time, eroded (typically by rainsplash) so that the material moved to the surface of the nest is then spread as a surface veneer over the soil surface, over an area that might be larger than that contained by the nest itself. The underground galleries of the nests are the material that is moved to the surface and may be quite extensive laterally (e.g. termite nests may extend 10’s of metres laterally) and vertically. Neither the lateral nor vertical density of the galleries is well defined (and likely to be species, climate and soil specific), nor is the spatial distribution of nests (all of which determine the mass of material moved). This also applies to the length of time before the galleries collapse and the soil above them subsides (which determines the rate of change of the porosity of the soil in those layers that contain the galleries).

The general principle here is that soil with a given size grading range (though the evidence for size selectivity for ants is not definitive, and nest-building behaviour appears to be strongly species dependent; Lobry de Bruyn and Conacher, Reference Lobry de Bruyn and Conacher1990, Reference Lobry de Bruyn and Conacher1994) is removed from each subsurface layer and deposited on the surface. To ensure that the mass in each layer does not change, we transfer material down from each overlying layer into the layer below until we reach the bottom of the bioturbed region. First, we will look at how to calculate the amount and grading of material transported to the surface by the ants. So for layer i the material removed is

$\Delta {\underset{¯}{g}}_i=B{\phi }_{i}A{\underset{¯}{g}}_i$(7.46)

where $\Delta {\underset{¯}{g}}_i$ is the change in the soil grading vector in this time step, B is the rate constant determining the total amount of material transported to the surface (not to be confused with B the matrix), φ_i${\phi }_i$ is a normalised function of depth describing the relative rate at which material is transported from each layer i (i.e. φ_i${\phi }_i$ is the depth dependency of ant-turbation) and matrix A is a normalised matrix describing the size selectivity of the ant-turbation mechanism. An example of matrix A is

$A=\left(\begin{array}{cccccc}{a}_{1,1}& 0& 0& 0& 0& \cdots \\ 0& a_{2,2}& 0& 0& 0& \cdots \\ 0& 0& a_{3,3}& 0& 0& \cdots \\ 0& 0& 0& 0& 0& \cdots \\ 0& 0& 0& 0& 0& \cdots \\ ⋮& ⋮& ⋮& ⋮& ⋮& \ddots \end{array}\right)$(7.47)

where the diagonal elements a_{j, j}$a_{j,j}$ are the proportion of the total mass in grading size range j that is transported upwards by the bioturbation (i.e. a_{j, j} = 0$a_{j,j}=0$ if nothing is removed, a_{j, j} = 1$a_{j,j}=1$ if it is all removed). In this example a_{4, 4} = 0$a_{4,4}=0$ because for grading size range 4 it is assumed that none of that size range is transported upward (e.g. because it is too coarse). I have shown only the first five rows and columns of the matrix for brevity.

It also turns out to be convenient to define the total mass (i.e. across all grading size classes) in layer i, m_i, and the total mass of the change in the soil grading vector for layer i, Δm_i,$\Delta m_i,$

$m_i=\sum _k{g}_{i,k}$(7.48)

$\Delta m_i=\sum _k\Delta g_{i,k}$(7.49)

where the subscript i indicates layer i and the subscript k is the kth grading class in layer i.

For simplicity of explanation let us consider a case where the ant galleries are excavated in the top four layers, and the layers below are untouched. Then material will be excavated from layers 2 through 4 and deposited into layer 1. To ensure that the mass in each of the layers is unchanged, we then need to transfer the excess material in top layer down into layer 2, from layer 2 to layer 3 and so on. (Figure 7.9). The mass and grading of the material brought to the surface is

$\sum _{i=2}^{4}B{\phi }_{i}A{\underset{¯}{g}}_i$(7.50)

The fluxes of soil between the layers are given in Figure 7.9 and are based on a fully explicit numerical representation of the process where each of the soil layers is fully mixed. Importantly in Figure 7.9 the grading of the material from one layer to the layer below is based on the grading of the layer prior to the redistribution of the soil into the layer above. This is the explicit numerical approximation of the time discretisation. There is no conceptual reason why the derivation can’t use the grading in the layer above after the soil has been redistributed to it; the arithmetic in the equation below is just more complex.

Figure 7.9: Schematic of how layers interact in ant and termite bioturbation. While the calculations are done using mass, the figure shows the layer thicknesses: (a) the soil profile before ant bioturbation, including the thin surface armour layer, (b) the movement of material by ants, from the top three layers into the surface armour layer (those layers above the thick line are subject to ant transport, those below not; the grey region is the surface armour layer; note that the only layer that changes thickness is the surface armour layer, so the bulk density in the lower layers is reduced), (c) the redistribution of excess sediment from the surface armour layer (the thickness of the lines and arrows indicate quantities of soil being moved between layers) and (d) the reinstatement of the armour layer, which is now a mix of the original armour layer and the ant transport material. The soil surface has risen because of the reduction of the bulk density of the layers from which ants removed materials.

We can write the equations for the new grading in each of the four layers. Starting at the surface and working down the profile,

$\begin{array}{c}{\underset{¯}{g}}_{1,\text{new}}={\underset{¯}{g}}_1+0-\left(\frac{\sum _{i=2}^4\Delta m_i}{m_1+\sum _{i=2}^4\Delta m_i}\right){\underset{¯}{g}}_1+\sum _{i=2}^{4}B{\phi }_{i}A{\underset{¯}{g}}_i\\ {\underset{¯}{g}}_{2,\text{new}}={\underset{¯}{g}}_2+\left(\frac{\sum _{i=2}^4\Delta m_i}{m_1+\sum _{i=2}^4\Delta m_i}\right){\underset{¯}{g}}_1-\left(\frac{\sum _{i=3}^4\Delta m_i}{m_2+\sum _{i=3}^4\Delta m_i}\right){\underset{¯}{g}}_2-B{\phi }_{2}A{\underset{¯}{g}}_2\\ {\underset{¯}{g}}_{3,\text{new}}={\underset{¯}{g}}_3+\left(\frac{\sum _{i=3}^4\Delta m_i}{m_2+\sum _{i=3}^4\Delta m_i}\right){\underset{¯}{g}}_2-\left(\frac{\Delta m_4}{m_3+\Delta m_4}\right){\underset{¯}{g}}_3-B{\phi }_{3}A{\underset{¯}{g}}_3\\ {\underset{¯}{g}}_{4,\text{new}}={\underset{¯}{g}}_4+\left(\frac{\Delta m_4}{m_3+\Delta m_4}\right){\underset{¯}{g}}_3-0-B{\phi }_{4}A{\underset{¯}{g}}_4\end{array}$(7.51)

where the first term on the right-hand side of each equation is the grading before ant-turbation, the second term is the soil being transferred down from the layer above, the third term is the soil being transferred down to the layer below and the fourth term is the material moved by the ants/termites. The zeros in two of the equations are placeholders so that the third and fourth terms are the same processes in all four equations. Constructing the transition matrix from these equations involves filling the matrix according to Equation (7.51):

$B=\left(\begin{array}{cccccc}-\left(\frac{\sum _{i=2}^4\Delta m_i}{m_1+\sum _{i=2}^4\Delta m_i}\right)I& B{\phi }_{2}A& B{\phi }_{3}A& B{\phi }_{4}A& \left(0\right)& \cdots \\ \left(\frac{\sum _{i=2}^4\Delta m_i}{m_1+\sum _{i=2}^4\Delta m_i}\right)I& -\left(\frac{\sum _{i=3}^4\Delta m_i}{m_2+\sum _{i=3}^4\Delta m_i}I+B{\phi }_{2}A\right)& \left(0\right)& \left(0\right)& \left(0\right)& \cdots \\ \left(0\right)& \left(\frac{\sum _{i=3}^4\Delta m_i}{m_2+\sum _{i=3}^4\Delta m_i}\right)I& -\left(\frac{\Delta m_4}{m_3+\Delta m_4}I+B{\phi }_{3}A\right)& \left(0\right)& \left(0\right)& \cdots \\ \left(0\right)& \left(0\right)& \left(\frac{\Delta m_4}{m_3+\Delta m_4}\right)I& -B{\phi }_{4}A& \left(0\right)& \cdots \\ \left(0\right)& \left(0\right)& \left(0\right)& \left(0\right)& \left(0\right)& \cdots \\ ⋮& ⋮& ⋮& ⋮& ⋮& \ddots \end{array}\right)$(7.52)

where, for brevity, only the first five rows and columns of the transition matrix are shown.