use the steady state and the non steady state approach respectively.
2.4.1 Steady state approach
This approach is used if the total inventory and the depth distribution of 7Be at the study sites and the reference sites have uniform distributions (steady state) before the erosion occurs. After an erosion event the inventory of 7Be (Bq m-2) at study sites is measured and compared to that at the reference site and eroded depth is calculated by following method:
If z is the layer of soil eroded at the study plot then the inventory of 7Be at the sampling point after erosion can be calculated by integrating (9) from z to .
(19) by rearranging (19) and using (11) we get,
§ ·
Δ = λ ⋅ ¨© ¸¹
Re f s
D A
z ln
A ȋʹͲȌ
If z in (20) is negative it implies that erosion has occurred and vice versa.
41
2.4.2 Non steady -state approach
This approach is used for calculating the eroded depth z for multiple erosion events. In situations where the erosion events are separated by a short time interval
Δt, the steady state approach is not applicable because of large differences between the reference inventory and the inventory of 7Be at the measurement point. This difference is due to the loss of 7Be by erosion at the sampling point after the first event and not enough time for the 7Be to reach a steady state (~106 days [33]). For such situations we cannot compare the inventories of 7Be at a study plot with that at the reference site.
2.4.2.1 Crank-Nicolson scheme
For non steady-state conditions the numerical Crank-Nicolson scheme was used to estimate the 7Be inventory within short time intervals between two erosion events.
This is a standard scheme to solve diffusion equation, as it is unconditionally stable and highly accurate. It uses a finite differences approach for the discretization of differential equations. Its application to the diffusion model given by equation (5) is given by equation C1. A detailed discussion on the Crank-Nicolson scheme is given in Appendix D.
− +
− − −
− +
§ λ·
− Δ +¨©Δ +Δ + ¸¹ − Δ
§ λ·
= Δ +¨©Δ −Δ − ¸¹ + Δ
m m m
i 1 i i 1
2 2 2
m 1 m 1 m 1
i 1 i i 1
2 2 2
D 1 D D
C C C
t 2
2 z z 2 z
D 1 D D
C C C
t 2
2 z z 2 z
(C1)
The MATLAB implementation of equation (C1) is given in Appendix F. To test the ability of this numerical scheme to calculate inventories of 7Be in soil the following simulations were performed.
1 Comparison of numerical and analytical solution
The Numerical solution was compared with the steady state solution and the truncation errors determined for the approximated solution. The transient behaviour of the solutions to equation (5) is smooth and bounded, meaning the solution does not develop local or global maxima that are outside the range of the initial data that is, the Crank-Nicolson scheme is unconditionally stable. The time increments (t)
42
and the depth increments (z) were chosen to obtain the smallest difference between the numerical and analytical solution. The percentage error for this combination of M and Nwas < 1%. Depth and time increments, z = 10-5 m and t= 0.06 day were estimated for tmax = 90 days and z= 2 cm. The simulations were performed for a number of time steps, M = 1000, 2000 and 2500 and a number of space steps, N = 500, 1000 and 2000 used for the solution matrix. For each combination of M and N, the numerical solution was compared with the analytical solution and absolute and percentage errors were determined. The best estimation for depth and time increments for the minimum error was N = 2000 and M= 1500.
2 Truncation error analysis for the Crank-Nicolson scheme
The truncation error for the Crank-Nicolson scheme is σ Δ( z )2 + σ Δ( t )2 ǡ
represents the rate at which the truncation error goes to zero.
Figure 12: A: Truncation errors for Crank-Nicolson scheme as a function of x for fixed t. B: Truncation errors for Crank-Nicolson scheme as a function of t for fixed x.
The true magnitude of the Truncation Error (TE) is given as:
= Δ +t 2 xΔ 2
TE K t K z
43
Kt and Kx are constants that depend on the accuracy of finite difference approximations.
To make TE arbitrarily small, both z and t must approach zero. Figure 12 gives the comparison of truncation error for Crank-Nicholson scheme for fixed time and space steps, t and x. For a fixed time step t, figure 12A shows that the dependence of a truncation error for Crank-Nicolson scheme is nearly identical except for small x. Figure 12B shows how the truncation error of the Crank-Nicolson codes varies when t is reduced and x is fixed. The solutions were obtained for diffusion coefficient, D = 3.65 × 10-13 m2 s-1, depth of the soil column L= 2 cm and tmax = 90 days. The Crank-Nicolson scheme can obtain solutions for much larger t as it is unconditionally stable. As Δ →t 0 for finite x, the TE of Crank Nicolson scheme approaches a constant value.
2.4.2.2 Erosion quantification using the Crank-Nicolson scheme
The erosion estimation technique for non steady state condition can be understood by considering two erosion events E1 (t1) and E2 (t2) occurring at time t=t1 and t=t2
(Figure 13).
Figure 13: Scheme for the calculation of erosion rates for multiple erosion events. E1 and E2 are erosion events at time instances t1 and t2 respectively. t’2 is the time instance before erosion event E2 occurs.
The necessary input parameters to solve the differential equation (5) using the numerical Crank-Nicolson method are as follows,
44
1. The initial condition is chosen as the inventory of 7Be after the erosion event at t=t1 C(z, t=0) =CE as per equation (19). CE is obtained from the measurement of
7Be inventory at a sampling point after an erosion event.
2. The input flux I0 (Bq m2 s-1) is estimated from the reference inventory by using equation (11).
3. The diffusion constant (m2 s-1) is evaluated by fitting equation (8) to the depth distribution of 7Be at the reference sites.
Using the three input parameters mentioned above, depth distribution profile of
7Be was numerically simulated until the time t =t'2. The simulated profile at t = t’2 of
7Be is integrated over the depth to obtain the total inventory (equation (21)). The simulated total inventory at point t = t’2 is denoted by ANum.
zEnd
'
Num Num i 2 i
0
A =
³
C (z , t=t )dz (21)here zi = − Δ(i 1) x, i 1, 2,...N ; N is the total number of spatial nodes including = those at the boundary and x is the spacing between zi. In equation (21), ANum is the total inventory of 7Be before erosion in Bq m-2 and CNum is the simulated volumetric concentration of 7Be in Bq m-3.
§ ·
Δ = λ ⋅ ¨© ¸¹
Num s
D A
z ln
A (22)
The inventory of 7Be obtained from (21) is compared to the inventory of 7Be measured at the study site after erosion event E2 at t = t2 and eroded/deposited layers are calculated by using (22).
45