- Email: [email protected]
Modelling phosphorus transport in soils and groundwater with two-consecutive reactions
Wat. Res. Vol. 25, No. 10, pp. 1205-1216, 1991 Printed in Great Britain.All rights reserved
0043-1354/91 $3.00+ 0.00 Copyright © 1991 PergamonPressplc
MODELLING PHOSPHORUS TRANSPORT IN SOILS A N D GROUNDWATER WITH TWO-CONSECUTIVE REACTIONS SUPRI NOTODARMOJO1, G. E. Ho ~O, W. D. SCOTT~and G. B. DAVIS2 tSchool of Biological and Environmental Sciences, Murdoch University, Murdoch, WA 6150 and 2CSIRO Division of Water Resources, Private Bag, P.O. Wembley, WA 6014, Australia (First received September 1990; accepted in revised form March 1991)
Abstract--A model is presented for one-dimensional transport of phosphorus (P) in soils and groundwater. Convective transport, hydrodynamic dispersion and time-dependent phosphorus sorption are accounted for in the model formulation. Time-dependent sorption of soil-P is considered to follow the two consecutive reaction model of Barrow and Shaw (J. Soil Sci. 30, 67-76, 1979) which has been extensively tested against experimental data and can be described by S = k • Cnt ~. The assumed sorption model allows parameters to be obtained by independent batch and column experiments. A numerical technique is used to solve the solute transport equation incorporating a correction to numerical dispersion to improve the numerical solution. An analytical solution for a simplified case is also presented to test the numerical technique. Parameter sensitivity analysis shows that influent concentration and the parameter k strongly affect the initial breakthrough time of solute, with m and n affecting the shape of the breakthrough curve. Preliminary investigations show that the applicability of the model to describe column experimental breakthrough curves is promising. Key words--phosphorus transport, time-dependent sorption, two-consecutive reactions, solute transport
modelling
INTRODUCTION
Understanding the movement of inorganic phosphates in soil is essential to predict the potential of phosphorus pollution of groundwater from fertilizer application. Management of the problem can be greatly assisted if a model is available that can quantitatively describe the movement of phosphate with infiltrating water and its interactions with soil. Phosphorus contamination of groundwater from application of inorganic and organic fertilizers has in many cases resulted in euthrophication of receiving waters when the groundwater discharges to shallow surface water bodies. The dangers of excessive phosphorus leaching are well documented (Birch, 1982; Lukatelich et al., 1987). Reactions of phosphates with soils are time-dependent (Barrow and Shaw, 1974; Munns and Fox, 1976; Ryden et al., 1977) and several models have been proposed. The reactions are generally viewed as either phosphate adsorption by or resorption from the soils. Selim et al. (1976a) combined nonlinear-kinetic and instantaneous terms, Cameron and Klute (1977) developed a model using a combination of linear-kinetic and linear-instantaneous terms. These models assume that soil particles consist of two different reaction sites. De Camargo et al. (1979) also used the two-site model. Their model was similar to Cameron and Klute's (1977) model, but they assumed that the kinetic-type reaction was due to adsorption sites which were less accessible. Mansell et al. (1985) modelled phosphate transport in a sandy soil by using
the three compartment model of Barrow and Shaw (1975). The model assumes that sorption-desoprtion follows two consecutive steps. During the first step, physical adsorption rapidly removes phosphate from solution, and during the second step the adsorbed phosphate undergoes chemisorption following an n th order kinetic reaction. Recently, Gerritse (1989) derived a phosphate transport model for sandy soils incorporating a time-dependent reaction. The reaction of phosphate and soil in his model was assumed to proceed according to the first-order rate equation of Munns and Fox (1976). The above models have been tested against experimental data and are potentially applicable to describe phosphate transport in soils. Generally, the models include kinetic terms with parameters which are difficult to obtain independently. In the models proposed by Selim et al. (1976a), Cameron and Klute (1977) and De Camargo et al. (1979), the proportion of the two different sites to produce a good fit of the model to the observed breakthrough experimental data has to be obtained by trial and error, while for the models of Mansel et al. (1985) and Gerritse (1989), the experimental breakthrough curve is needed to determine the ratio of forward and backward physical adsorption. It would be desirable to develop a solute transport model which allows an independent evaluation of the sorption-desorption parameters. Although there have been many attempts to model the time-dependent sorption of phosphate by soils, the model proposed by Barrow and Shaw (1979) appears to have a theoretical basis and fits
1205
1206
SuPra NOTODARMOJOet al.
experimental data of phosphate sorption by soils from many parts of the world (Probert, 1983; Barrow, 1989). No transport model appears to have been developed incorporating this more universally applicable sorption equation. The objective of this paper is to simulate the movement of phosphate in soils and groundwater in one-dimensional steady flow conditions, incorporating concentration and time dependent reactions, and with model parameters that can be evaluated independently from batch and column experiments. THEORETICAL
The two-consecutive reaction model o f P sorptiondes.orption
Sorption experiments are often conducted over a short period of time and equilibrium is assumed between sorbed and liquid phases. For phosphorus, true equilibrium is not reached within a short period. The reaction of phosphate with soils is complicated and the process is still puzzling. Much research has been undertaken to elucidate the processes involved. A plausible and likely mechanism has been suggested by Barrow and Shaw (1975) consisting of a threecompartment model for phosphate reaction with soils: A ~--,B ~--~C. Compartment A contains phosphate in the soil solution. Its concentration is determined by the concentration of adsorbed phosphate in compartment B according to an adsorption process. Compartment C contains phosphate which has been converted into a form which is not in direct equilibrium with compartment A. The reaction of dissolved phosphate (in compartment A) with the soil is considered as a rapid surface adsorption (or desorption) onto reactive surfaces followed by a slow penetration (or release) by solid state diffusion of this surface phosphate into subsurface horizons within the interior of the particles (compartment C). Barrow and Shaw (1979) concluded that the best description of the combined effect of both on phosphate sorption is of the form: S =k.C".t m
(1)
where S is the phosphate adsorbed onto the soil particles (~g/g of soil), C is solute concentration in the pore space (#g/cm3), k is a constant, which gives a relative measure of the adsorbent capacity of the soil for phosphate, n and m are constants which indicate the dependence of sorption on concentration and time (t) with a high value indicating stronger dependence. Barrow (1983) showed this simple descriptive model to be in good agreement with a more complex mechanistic model. A plot of log C vs log S from a series of batch experiments conducted at various concentrations and times of incubation would give the constants k, m and
n. The constants are affected by many factors such as: properties of the soils, pH, temperature, electrolyte background and soil-to-solution ratio (Barrow and Shaw, 1979; Syers and Iskandar, 1981). Phosphorus transport with two-consecutive reactions
In many cases the horizontal movement of pollutants such as phosphate in soils and groundwater can be simulated by a one-dimensional model (Selim et al., 1976b; Kinzelbach, 1986): [ O C ] ~ q C O S OOCst =--OzO Ov.O .~z - - ~ z - p "-~
(2)
where q is the volumetric flux of water through unit cross sectional area (cm/day), v = q/O is the average pore water velocity (cm/day) and Dv is the hydrodynamic dispersion coefficient (crn2/day). S is the solute associated with the solid phase (#g/g of solid), p is the bulk density of solid, C is the solute concentration in the solution phase (ktg/cm ~) and 0 is the volumetric water content. For steady uniform water movement the variables q, 0, v and D v are constant. Equation (2) then can be simplified to: OC 8t
OzC ~C D~. sz 2 - v c3z
p ~3S 0 Ot
(3)
Assuming S can be described by equation (1), then: 8S k . C , m.t,,_~ + k . t , , . n . C , _ l O C " dt = Ot
(4)
Substitution of equation (4) into equation (3) results in: OC c~2C - - =Dv.~z2 - V ~t
OC
"~ --~P . k . C " . m . t -o.[k'tm.n'C"-
m-l
i
8C
1~-.
(5)
Letting R = l + K ' t " ' n - C n ] and K = p / 0 " k , and rearranging equation (5) we then get OC
Dv 02C
O---t= R
dz z
v ~C
K
R Oz
R C"'m'tm-I"
(6)
Equation (6) is the solute flow equation with time-dependent sorption assuming no stagnant (immobile) water phase within the pores of the media. Analytical solution for a simplified case
Because of its non-linear nature, equation (6) is difficult to solve analytically. However, simplifying the equation allows analytical solution. When the value of parameters rn and n are equal to 1, equation (4) becomes: ~S c3t
~C 0t
-- =k.C +k.t.--
(7)
and equation (6) becomes: OC 02C 0C (1 + K . t ) - ~ = O~ . -~-fz2 -- v .-~z - K . C
(8)
Phosphorus transport in soils and groundwater Assuming a semi-infinite column length, a well mixed feed solution at the top boundary (z = O) and initially no sorbed or liquid phase phosphorus in the column, the appropriate boundary and initial conditions become: C=Co
at
z--O
~C dz
----0
as
C=S=0
at
t--0
for allz.
dC O2C OC K ' z "-~z = D "-~z 2 - V "-~z
(9)
-
K "C
(10)
with transformed initial condition: T=I.
(11)
Further using the Euler transform for z, z=C,
then
dz=e".dt?
(12)
and setting C(z, z) = ok(z, rl)
(13)
gives: K . ~ n~
C(z,t)
- 02~b 0~b = uv.-~z: -- V.~z - K . q b
at
t?=0;
a¢
--~0 Oz
q~=C0
as
f F R.z-v.t
z--,~.
at
(14)
z=0
Co f
(17)
Numerical solution for general cases When analytical solution is not possible, numerical approximation can often be used successfully. For one-dimensional solute transport, finite difference techniques are widely used. The explicit finite difference technique is commonly used, even though it may require extended computing time because of its restrictive stability criteria. However, by introducing corrections for numerical dispersion (Appendix A), the stability criteria become less restrictive. The (z, t) plane is discretized by a rectangular grid of points: i = 1, 2, 3 . . . . , I along the z-axis, and j = 0, 1, 2, 3 . . . along the t-axis.
(15)
The distance between the nodes Az remains constant:
Equation (14) together with the initial and boundary conditions [equation (15)] has been solved by van Genuehten (1981) to give the analytical solution: C(z,t)=T'~exP~
q
Figure 1 illustrates concentration breakthrough curves calculated by using equation (16) for various values of K. The number of pore volumes is defined as V/Vo where V is the total volume of liquid phase which has flowed through the column at a given time and V0 is the volume of liquid phase in the column. It can be observed that the curves exhibit a fairly flat top because of the strong dependence of sorption on time. A high value of parameter K gives a low final quasi-steady outlet concentration and a much delayed breakthrough.
with the following initial and boundary conditions: q~=0
i
Co = .er
2-7(B.-~Tt-)~,2 ].
z--.oo
at
column is available (Lapidus and Amundson, 1952):
for a l l t
Now if we let z = 1 + K . t , then d z / K = d t and we can rewrite (8) as:
C=0
1207
/(v
-
Az = L / ( I - 1)
~)z'x
2-~.-~ )
100--
x erfc; (K.z - ,,,,) \ 2 ' ~ 1
dT--
(18)
+ e x p/(v t--t \ 2-D, ,/
/j
x erfc( ( K . z + #~1) ~
tO
80-
0 tO X O3
(16)
K=0.2 K=0.4 K = 0.6K = 0.8
,f K=I K=2
D v : 25 cm2/d 1O0 cm/d
60--
t13 t~ p.
where #=v.
.o
1
+ 4 " K ' D , ~ 1/2 v----T--- ]
40--
cO3 °E -D
20--
rl = ln(l + K. t) and effc(x) is the complementary error function. With the sorption-desorption simplified to follow a linear isotherm, S = k . C, an analytical solution to the solute transport equation for a concentrationtype inlet boundary condition and semi-infinite
0.0
I 0.5
1.0
1.5
2.0
2.5
3.0
N u m b e r of pore volumes
Fig. I. The effect of parameter K on the breakthrough curves calculated with equation (16).
SUPRi NOTODARMOJOet al.
1208
where L is the column length and I is the total number of the nodes along the column. The time difference At is determined by the stability criteria of the scheme. The finite difference equation is derived by using a Taylor expansion (Mickley et al., 1957; Carnahan et al., 1969). Here we use a forward difference in time, a central difference for the second derivative, the dispersive term, and a backward difference in space for the advective term: 0C ~t
Ci +l - C ! At
Dv [C~+a - 2"C~ + CJ-I? R Az 2 e ~C~-C~:_,~ K t ~(19) R" L _ R" [Cj],.m . 1 where
R = l + K'tm'n'[C~]"
~.
(20)
Equations (19 and (20) together with appropriate boundary and initial conditions [equation (9)] can be used for the computation; however, a correction to the hydrodynamic dispersion is necessary to minimize the inherent differencing error and to improve the performance of the finite difference technique. Van Genuchten and Parker (1984) showed that with the concentration fixed at the inlet and a linear adsorption isotherm, a better result was gained if the
column was considered as a portion of an infinite column. The exit concentration is then simply the concentration at position z = L Figure 2 is a schematic diagram of the column system and the finite difference compartments.
Correction due to numerical dispersion Numerical approximation of a partial differential equation deviates from the exact solution. The deviation is a direct result of replacing the original differential equation [e.g. equation (6)] by a difference equation with a finite time step (At) and layers of finite thickness (Az). This deviation can be quite substantial and is called numerical dispersion (van Genuchten and Wierenga, 1974). It was first quantified by Chaudari (1971) as resulting from the neglect of higher order terms in At or Az in the Taylor series approximation, van Genuchten and Wierenga (1974) and De Smedt et al. (1981) incorporated more terms, and derived techniques to reduce the errors caused by numerical dispersion. They showed how these terms could be found in the case of classical hydrodynamic dispersion, in both the case of steady water flow and the case of transient unsaturated water flow with the presence of immobile water. Only the important terms of order At (the lower orders) were taken into account. By using a similar technique, a correction to the numerical dispersion is derived in our non-linear case (Appendix A) and the result can be written in the following form:
D* = Dv - ~ .(Az--AtR'V)
(21)
Co (influent concentration)
where D* is an effective hydrodynamic dispersion coefficient. When used in place of Dr, it should minimize numerical dispersion. 1 Top
of the c o l u m n
,=~---- 2 ,=-- 3
Az = Depth increment
Stability of finite difference scheme One method to test the stability of a finite difference technique uses a matrix method (Smith, 1978); the technique was used by Binh (1981). Using this method, the difference equation for the partial difference equation [equation (19)] is applied to all compartments in the system. Applying the Gerschgorin circle theorem (Appendix B) to this matrix, we find the stability criteria for the time-step is: At ~<
,,=,-- L-1
L
Bottom of the c o l u m n , Effluent concentrations L+I are computed at this point. L+2 Bottom of imaginary column, assumed at infinity, where concentratlon gradient is zero C e (effluent concentration)
Fig. 2. Representation of the column for numerical solution.
1
2.Dv v K R.Az 2 t - ~ . A z + ~ . R C n ' m ' t " - I
(22)
Another constraint that has to be satisfied is the mass balance of the system. The mass balance criteria can be expressed by the Courant and Neumann criteria (Kinzelbach, 1986). The Courant criterion requires that: Az At ~<--. v
(23)
In physical terms, the criterion states that no more
Phosphorus transport in soils and groundwater pollutant mass can leave a segment via convection during the time interval [t + At] than is inside at the beginning of the time interval. The Neumann criterion, which is analogous to the Courant criterion, ensures that the concentration gradients cannot be reversed by dispersion-diffusion fluxes alone. It has the form: ~LZ 2
At
(24)
~ < -
2-D*"
AND DISCUSSION
Validation of numerical solution The finite difference method is said to converge if the discretization error approaches zero as the grid spacings Az and At tend to zero (Carnahan et al., 1969). When the stability criteria are satisfied we have found that our numerical solutions are convergent. To assess the accuracy of the numerical solution, we compare numerical results to those from analytical solutions for the simple cases. This approach was used by van Genuchten and Wierenga (1976) in solving equations for the convective--dispersive trans-
100 - -
v ,. L ,n ,, m -
D
port of reactive solute with exchange between mobile and stagnant water. Figure 3 shows a comparison between the numerical approximation and the analytical solution as given by equation (16), while Fig. 4 shows a comparison with the analytical solution derived by Lapidus and Amundson (1952). It can be observed that the agreement between numerical approximation and analytical solution is very good in both cases.
Model sensitivity analysis
For a given calculation cycle, the value of At obtained with equations (22) and (24) will vary for every compartment. The minimum calculated At for all compartments is then used as the criterion for the next step. For most cases, equation (22) has been found to be adequate to ensure stability. However, for solute transport where the convective flux is dominant or where the solute exhibits a very strong reaction with the media, a combination of equations (22), (23) and (24) needs to be used. RESULTS
1209
100 c m / d 2 0 cm 1 1
To evaluate the sensitivity of the solution of the phosphorus transport model to different model parameters, breakthrough curves are calculated for a 20 cm long soil column leached with 3 pore volumes of phosphate solution followed by solute free water at a constant flow-rate. The results of this simulation are presented in Figs 5-9. The influent concentration has a profound effect on the effluent concentration (Fig. 5). The retardation or "hold up" time is large for low influent concentrations, and for high concentrations, the effect of sorption is reduced and the shape of the c u r v e converges to that of a non-reactive solute case, with an S-shaped curve. When the incoming solute concentration is low, the gradient for transfer into the grains becomes low, hence the process of transfer to the solid occurs only at a very low rate. This does not lower the concentration at the surface and consequently lessens the surface reaction rate. This mechanism explains why at lower input concentrations applied to the soil, the desorption curve shows more tailing. When K increases, the hold up time also increases significantly although not in a linear proportion, and
D v" 5 cm2/d K - 0.5
NA
80-0
C 0
.--_ x u:,
--
:
60
NB
--
D v" 50 cm2/d K-2.5
u) (I)
-~
NA
0
"~ C
D v - 100cm2/d
4 0 --
K ,, 5 . 0
O'-E - NB~
20-
0
1
I
I
I
I
I
2
3
4
5
6
Number
of pore volumes
Fig. 3. A comparison between analytical solution and numerical solutions. AN = analytical solution, NA = numerical solution with correction and NB = numerical solution without correction.
1210
SUPRI NOTODARMOJOet al. 100
cm
~~'~=-~2
v = 100 cm/d
--
L = 20 n =l m-0
NA
v
80--
O tO O
N
/d
K=0.5
X
60--
Dv=5Ocm2/d K=1.5
t-
.2 tO
40--
E
~5 20--
NB 0
1
2
3
4
N u m b e r of p o r e v o l u m e s
Fig. 4. A comparison between analytical for n = 1 and m = 0 and numerical solutions. LA = analytical solution, NA = numerical solution with correction and NB = numerical solution without correction. the peak concentration decreases concurrently (Fig. 6). This is a direct consequence of the parameter K being a relative measure of the adsorbent capacity for phosphate and proportional to the reactive surface of the soils where the instantaneous reaction between soil particles and solute occurs. The magnitude of parameter n for the soilphosphate interaction is expected to be less than unity. The parameter n has a significant effect on the starting time of the breakthrough curves for low n values (Fig. 7). However, the slope of the
breakthrough and the descending portion of the curves are most affected by parameter n. The value of n relates to the availability of surface reaction sites and governs the ability of the particles to adsorb solute at higher concentrations. The slope of the curve is proportional to the ratio of parameter n to m. When n/m is close to unity, the shape of the curve is steeper since the rate of solid diffusion, which is responsible for the time dependency, is equal to the surface reaction rate at a given concentration.
Dv=25cm2/d v = 50 cm/d K=1.2 n=0.7 m=
100
8O O tO O
c 0
60
--
40
--
E
~5 20 --
0
1
I
I
I
,
2
3
4
5
Number of pore volumes Fig. 5. The effect of influent concentration on the breakthrough curves.
6
Phosphorus transport in soils and groundwater
1211
K =0.2
100
K=0.4 K =0.6 A
8O
K =0.8
o~ v
K=I.0
O ¢O O x (D tO tO O
.£
K=2.0
60
K=4.0
40
tO ri f
E
~5 20
I 0
I 1
2
3
4
5
6
Number of pore volumes Fig. 6. The effect of parameter K on the breakthrough curves. An increase in m lowers the peak of the curves (Fig. 8). The parameter m affects the potential of phosphate storage in the interior of soil particles; a high value increases this potential and means that more can be absorbed during the sorption period and more can be leached from the soil particles during the leaching process. The effect of dispersion is illustrated in Fig. 9. Dispersion affects the starting point of the breakthrough, with a relatively little effect on the shape of the curves.
~
100 - -
A
80--
o~ v 0 r0 0 .~_ x q~
The general pattern of the relative exit concentration distribution is a non-symmetrical curve, with a sharp front in the ascending portion, followed immediately by a decrease in effluent concentration after displacement with solute-free solution, then a significant tailing takes place. The net effect of a change in one parameter is not independent of the values of other parameters. The extent of tailing has important implications for the phosphate decontamination of soils and groundwater. Once contaminated by excessive phosphates, a
60--
n=l.0 n=0.9 n=0.8 n=0.7 n=0.6 n = 0.5
to
C
40
C
E ~5
20
0
1
2
I 3
4
N u m b e r of p o r e v o l u m e s
Fig. 7. The effect of parameter n on the breakthrough curves.
5
6
1212
SUPRI
100 -
NOTOD.CtMOJOet al.
D v = 25 cm 2/d v = 50 cm/d K=1.2 C O = 10 mg/I n=0.7
80 A
o~ 0 0 t.)
60--
x o') ~D
._(2
40--
r,
E ~5
20--
1
2
3
4
5
6
N u m b e r of p o r e v o l u m e s
Fig. 8. The effect of parameter m on the breakthrough curves. soil may take considerable time to desorb or displace sorbed phosphates due to the time-dependent sorption characteristics of phosphate in soil. Even with relatively low concentrations, the sensitivity analysis shows that the displacement of phosphate from soils may occur over an extended period. An example o f model application To illustrate the practical application of the model, two breakthrough curves from two different soil types are presented. Batch experiments were carried out to obtain estimates of the parameters k, n and m, according to Barrow and Shaw's (1979) method. Leaching experiments were carried out by using 20 cm soil columns following the method of Mansell et al. (1985). Hydrodynamic dispersion coefficients v = 100 cm/d K-1.2 C o = 10 mg/I n=0.7 rn = 0.4
100 -
80-
were also determined by using a similar miscible displacement experiment technique. Details of experiments and more extensive results will be presented in a separate paper. Figure l0 shows the breakthrough curve from an Australian sandy soil (Gavin A horizon), for which the important properties were given by McArthur and Bettenay (1960). Figure 11 shows the breakthough curve from an Indonesian loamy clay soil (Buah Batu soil). As can be observed, the agreement between experimental results and simulated curves are good. Further development We have developed a one-dimensional model for the transport of soluble phosphate incorporating time dependent sorption using a relationship that has been
~
Dv = 5 cmS/d
V. oc..
'
O tO 0
\\~--- D,=IO
60-
cm2/d
X (/)
C 0
4O
• Dv - 50 cm2/d
C
E
a
• Dv = 60 cma/d
2O
Dv
70
[
I
3
4
Cm 2 /dl
m o
2
~~
I 5
N u m b e r of p o r e v o l u m e s
Fig. 9. The effect of hydrodynamic dispersion on the breakthrough curves.
6
Phosphorus transport in soils and groundwater
A
100
--
80
--
1213
[ m
O e-" 0 O
D v " 15.50 cm2/d v - 6 3 . 1 5 cm/d Co - 2 . 5 0 mg/I K - 0.46 n = 0.46 m ,= 0 . 2 3
n
60
--
X
O
._o
40
--
20
--
t,--
E
a
O
cP []
I 0
1
I
I
I
I
2
3
4
5
Number
of p o r e v o l u m e s
Fig. 10. Experiment and simulated breakthrough curves from leaching with 2.50 mg/! phosphate solution for Gavin A soil.
tested over a wide range of soils and conditions. The model allows the independent determination of the sorption parameters. This is an improvement over similar existing models. A one-dimensional model can be applied in a wide range of field situations under restricted conditions where one-dimensional flow can be approximated. 100
Real systems are, however, generally very complex. The presence of immobile water cannot be ignored and a number of models referred to in the Introduction account for the time-dependent behaviour of sorption simply by its presence. Time-dependent sorption of phosphorus in soils occurs, however, even in batch systems. The presence of immobile water
m
cb
/° J
B
80
-
l
D v - 118.60 cm2/d v - 6 7 . 0 2 cm/d C O = 4 0 0 rng/I K= 197 n = 0.26 m=0.2
[]
oo
A v O cO O
m
60
-
40
--
X
._m C
O O
r0
t,-
O [] 0
E
a
20--
I 0
I
I
I
I
I
6
12
18
24
30
36
Number of pore volumes Fig. 11. Experiment and simulated breakthrough curves from leaching with 400 mg/l phosphate solution
for Buahhatu B soil.
1214
SUPRI NOTODARMOJOet al.
would add to this in a flow system, although this effect would be small when pore velocity is low. One-dimensional flow in the vadoze zone is characterized by varying soil water content. We are currently developing our model to incorporate the unsteady state condition due to irrigation and subsequent moisture redistribution on phosphorus transport. In particular we are investigating whether the steady state model can approximate the behaviour of the more complicated unsteady state model. Results of this study will be published separately.
CONCLUSION A non-linear phosphate transport model with two-consecutive reactions has been presented. The model can use independent batch experiments to obtain the value of the parameters that characterize the phosphate reaction with soils. An explicit finite difference code has been developed incorporating a correction for numerical dispersion. A new analytical solution for simplified conditions was derived and evaluated to aid validation of the finite difference code. Numerical results showed little numerical dispersion and were in close agreement with the new analytical solution and the analytical solution of a simple linear adsorption case. Parameter sensitivity analysis indicates that the model is very sensitive to the influent concentration (Co), and the parameter k, affecting the starting time of the breakthrough curve. Parameters n and Dv affect the start of the breakthrough curves although the effect of Dv is not as significant. The model also shows that at low applied phosphate concentrations, breakthrough may still occur, and sorbed phosphate can be desorbed by leaching. Further testing of the model against experimental data is planned Acknowledgement--We would like to thank International
Development Program of Australian Universities and Colleges (IDP) for providing a postgraduate scholarship to one of us (S.N.)
REFERENCES
Barrow N. J. (1983) A mechanistic model for describing the sorption and desorption of phosphate by soil. J. Soil Sci. 33, 733-750. Barrow N. J. (1989) The reaction of plant nutrients and pollutants with soil. Aust. J. Soil Res. 27, 475-492. Barrow N. J. and Shaw T. C. (1975) The slow reactions between soil and anions: 2. Effect of time and temperature on the decrease in phosphate concentration in the soil solution. Soil Sci. 119, 167-177. Barrow N. J. and Shaw T. C. (1979) Effects of solution and vigour of shaking on the rate of phosphate adsorption by soil. J. Soil Sci. 30, 67-76. Bear J. (1979) Hydraulics of Groundwater, pp. 225-276. McGraw-Hill, New York. Binh A. (1981) Neutralization of acid wastewater by limestone-bearing sand. Ph.D. thesis, Murdoch University, Western Australia.
Birch P. B. (1982) Phosphorus export from coastal plain drainage to the Peel-Harvey Estuarine system, Western Australia. Aust. J. Mar. Freshwat. Res. 33, 23-32. Cameron D. R. and Klute A. (1977) Convective~lispersive solute transport with a combined equilibrium and kinetic adsorption model. Wat. Resour. Res. 13, 183-188. Carnahan B., Luther H. A., and Wilkes J. O. (1969) Applied Numerical Methods, pp. 429-530. Wiley, New York. Chaudari N. M. (1971) An improved numerical technique for solving multidimensional miscible displacement equation. Soc. Petrol. Engng J. 11, 277-284. De Camargo A. O., Biggar J. W. and Nielsen D. R. (1979) Transport of inorganic phosphorus in an Alfisol. Soil Sci. Am. J. 43, 884-890. De Smedt F., Wierenga P. J. and van der Beken A. (1981) Theoretical and Experimental Studies of Solute Movement through Porous Media with Mobile and Immobile Water,
pp. 207-210. Vrije Universiteit Brussel, Belgium. van Genuchten M. Th. (1981) Non-equilibrium transport parameters from miscible displacement experiments. Research Report 119, USDA-ARS, Riverside, Calif. van Genuchten M. Th. and Parker J. C. (1984) Boundary conditions for displacement experiments through short laboratory soil columns. Soil Soc. Am. J. 48, 703-708. van Genuchten M. Th. and Wierenga P. J. (1974) Simulation of one-dimensional solute transfer in porous media. New Mexico Agric. Exp. Sm Bull. 628. van Genuchten M. Th. and Wierenga P. J. (1976) Numerical solution of the convective-dispersion with intra-aggregate diffusion and non-linear adsorption. In System Simulation in Water Resources (Edited by Vansteenkiste G. C.), pp. 275-292. North-Holland, New York. Gerritse R. G. (1989) Simulation of phosphate leaching in sandy soils. Aust. J. Soil Res. 27, 55~6. Kinzelbach W. (1986) Groundwater Modelling: an Introduction with Sample Programs in Basic, pp. 211-265. Elsevier, Amsterdam. Kreyzig E. (1983) Advance Engineering Mathematics, pp. 822-824. Wiley, New York. Lapidus L. and Amundson D. R. (1952) Mathematics of adsorption in beds. VI. The effect of longitudinal diffusion in ion exchange and chromatographic columns. J. Phys. Chem. 56, 984--988. Lukatelich R. J., Schofield N. J. and McComb A. J. (1987) Nutrient loading and macrophyte growth in Wilson Inlet, a bar-built south-west Australian estuary. Estuar. Coast. ShelfSci. 24, 141-165. Mansell R. S., McKenna P. J., Flaig E. and Hall M. (1985) Phosphate movement in columns of sandy soil from a wastewater-irrigated site. Soil Sci. 140, 594i8. McArthur W. M. and Bettenay E. (1960) The development and distribution of the soils of the Swan Coastal Plain, Western Australia. CSIRO Aust., Soil Publ. No. 16. Mickley H. S., Sherwood T. K. and Reed C. E. (1957) Applied Mathematics in Chemical Engineering. McGrawHill, New York. Munns D. N. and Fox R. L. (1976) The slow reaction which continues after phosphate adsorption: kinetic and equilibrium in some tropical soils. Soil Sci. Soc. Am. J. 40, 46-51. Probert M. E. (1983) The sorption of phosphate by soils. In Soils: an Australian View Point, pp. 427-435. CSIROAcademic Press, Sydney. Ryden J. C., McLaughlin J. R. and Syers J. K. (1977) Time-dependent sorption of phosphate by soils and hydrous ferric oxides J. Soil Sci. 28, 585-595. Selim H. M., Davidson J. M. and Mansell R. S. (1976a) Evaluation of a two-site adsorption-desorption model for describing transport in soil. Summer Computer Simulation Conference, 12-14 July 1976, Washington, D.C. Selim H. M., Mansell R. S. and Elzeftawy A. (1976b) Distribution of 2,4-D and water in soil during infiltration and redistribution. Soil Sci. 121, 176-183.
Phosphorus transport in soils and groundwater
1215
Smith G. D. (1978) Numerical Solution of Partial Differential Equations: Finite Differences Methods, 2nd edition. Oxford University Press, Oxford. Syers J. K. and Iskandar I. K. (1981) Soil-phosphorus chemistry. In Modeling Wastewater Renovation (Edited by Iskandar I. K.), pp. 575-599. Wiley, New York.
Further simplification by neglecting third order terms and non-linear terms in ac/at results in:
APPENDIX A
A comparison of equation (A7) with equation (6) suggests that an approximation to the effect of numerical dispersion from the forward time difference can be obtained when the difference equation uses an apparent or corrected dispersion coefficient: At .v 2 D*=Dv-~ 2.RZ. (A8)
Correction for Numerical Dispersion Consider the following Taylor series expansion of C{ + 1:
-C{+OC.Atm Cij+l _ Ot
02Catz 03CAt 3 + ~ T - - f -+ Ot3 6
+ terms with higher orders of At.
(AI)
Substitution of equation (6) into equation (A1) and rearrangement of the terms results in:
At
At
-
Oz2
R Oz
AC At - -
K =
----
v OC . C n . m . t m - l
R 3z
_[~
Oz2
O C
R Cn'm "tm-I
A C
Oz
J J C~+, Ci_, CiJ --
Az
2.Az
At 2 02 rDv 02C
v OC
K C".m.t m-I]
R az
R
AC At (A2)
R
v.Orl-]Oc
ot.Oz2
ot K
v O2C ROtOz
,
O [ ~ ' C "m'tm-I].
C{+l_C{
At
az 2
(A9)
-
1 '
~ t --K'[C~]"'m't"-t
_ vIc{-~c{_,
The second term on the right-hand side of equation (A2) can be considered as a first estimate of the effect of numerical dispersion (Chaudari, 1971; van Genuchten and Wierenga, 1974; De Smedt et al., 1981). Note that this only considers the effect of the approximation to numerical dispersion in time. This term is:
Dv 03C
J Ci_,
Az
+2
+ terms with higher order.
OF1 ] o 2 c
--
Az -C~+, - 2C~ + C~_ 1]
R Oz
+-6-~?[~ 0 :
At v2]02C
+T~J~z~ (A7)
It is known that central differences best approximate the spatial derivatives (van Genuchten and Wierenga, 1974). However, with backward differences in the advective term, equation (AS) has to be modified. Here we approximate the derivative with central differences and look for the error
AtO [ ~ 02C vOC K . c . . m . t m _ , ] + 2fit
- - - - - -
R
x
[C~+ t - 2C~ +
(AI0)
Comparing equations (6) and (AI0) the numerical model should use the following corrected hydrodynamic dispersion coefficient: (A3)
dt
van Genuchten and Wierenga (1974) and De Smedt et al. (1981) proceed as follows. As a first order approximation, only the advective term in equation (A3) is retained.
APPENDIX B
Stability Criteria of the Finite Difference Equation _ _
~ _ _
Ot
C n . m .l m-I
=
At
Oz2
R Oz
R
At O F
+-~
v OC]
L-~zj.
(A4)
Rearranging terms and expanding the last term, gives:
AC At
-
-
=
K C,.m.tm_ t R
--__.
v dC R Oz
At v 0 0C 2 R Oz Ot
DvO2C R c3z2 At OC O F l-]
AC -
=
At
--__
K.c,.m.tm-I R
Uj = A • U j - ]
(A5)
Considering equation (6) for OClOt and neglecting the last term of equation (A5): -
Finite-difference approximations contain errors which occur due to neglecting terms of higher order. These can often give rise to numerical inaccuracy or regions of solution instability. The stability of a finite difference scheme can be assessed using a matrix method proposed by Smith (1978). The finite difference representation of the governing partial differential equation (using backward differences for the advective term), with the top and bottom boundary conditions, is given in matrix form as:
v OC DvO2C R Oz2
(B1)
where U s is the vector of all values of Us at time j, U j is the vector at time j - 1, and A is the matrix containing all constants. By recursion technique, equation (B1) can be expressed as:
R 0z
UJ = A .U j-~ = A .(.4 . U j-2)
At v D~33C At v2 O2C Atv2 0C O [ 1 ] 2 R R Oz3 ~ 2 R 20z 2 + 2 R OzOz
(A6)
....
A ' A . . I A ( A . U °) (B2)
where U ° is the vector of all values at time j = 0. Then equation (B2) becomes:
U s = (A) j. U°.
(B3)
1216
S U P R I NOTODARMOJO et
which has the general form:
Defining the error vector e = U - U*, where U* is the vector of exact values, then: ~. = U j -
U * = ( A ) J ' ( U O - U * ) = AJ'Eo
C j÷ ~= A • CL
I/Jl ~< 1, k = 1,2,3 . . . . . n
Application of Gerschgorin's circle theorem to matrix A gives:
lit - ( 1 - 2 D - V - K ) I < ( D ) + ( D
Az
j
(BS)
(I-2D-
V-K)-(2D+
(1 - 4 D - 2 V - - K ) ~ #
D = ~-~7k'
v-
v'At
az.R'
<~1 - K .
(BII)
For stability: - I ~ - 1 only when: (1 - 4 D - 2 V - K ) > ~ I
(B6) or
Defining the following terms:
Dr'At
V)
or
J I
(B10)
~
C{_,l
R . [ C ]i , - , C ~.m.t,~
+ V)
where (1 - 2 D - V - K) are the diagonal elements (A,.,), and (D) + (D + V) is the sum o f moduli (2,) of the element along the row excluding the main diagonal for segments other than the boundary. Therefore,
where/~ is the eigenvalue of A (n x n). Now the problem is confined to testing the eigenvalues of matrix A. One method of finding the bounds of eigenvalues is Greschgorin's circle theorem (Smith, 1978; Kreyzig, 1983). This theorem states that if 2, is the sum of the moduli of the element along the s row excluding the diagonal element As.s, then each eigenvalue of A lies inside or on the boundary of a least one of the circles IV - As.~l = 2~. The finite difference equation (19) can be rewritten as:
R [_
(B9)
(]34)
where E0 is the vector of round-off error from the beginning. According to Smith (1978), a finite difference scheme is stable only when ~, remains bounded as j increases indefinitely. For this to happen the eigenvalues of A must have modulus values less than or equal to unity:
c~C _ C,J+' - C { _ D v f C { + , - 2 " C ~ + Ot At R L Az2
al.
4D + 2V + K ~< 2.
K=K'At[c~],_].m.t,,_l
Substituting the values of D, K and K into equation (B12) results in:
R
equation (B6) becomes: At ~<
c ? ' = [D + v]. c~_~ +[1-2.D-V-K].C~+DCJ,+~.
(BI2)
1
2.D v
v
(B13)
K
RAz E"FR~-'~z +--2"R "Cn-I "m "t m-I
(B7)
In matrix form this is: CI
c~ c~
)+1
1
CI
(D+V)
(1-2D-
( D + V)
V-K)
(I-2D-
c:
D
V-K)
D
(B8) w
m
c~ . Cn+
( D + V)
(I-2D-
( D + V)
V-K)
D
(1-2D-V-K)
c~ _Cn+l
0043-1354/91 $3.00+ 0.00 Copyright © 1991 PergamonPressplc
MODELLING PHOSPHORUS TRANSPORT IN SOILS A N D GROUNDWATER WITH TWO-CONSECUTIVE REACTIONS SUPRI NOTODARMOJO1, G. E. Ho ~O, W. D. SCOTT~and G. B. DAVIS2 tSchool of Biological and Environmental Sciences, Murdoch University, Murdoch, WA 6150 and 2CSIRO Division of Water Resources, Private Bag, P.O. Wembley, WA 6014, Australia (First received September 1990; accepted in revised form March 1991)
Abstract--A model is presented for one-dimensional transport of phosphorus (P) in soils and groundwater. Convective transport, hydrodynamic dispersion and time-dependent phosphorus sorption are accounted for in the model formulation. Time-dependent sorption of soil-P is considered to follow the two consecutive reaction model of Barrow and Shaw (J. Soil Sci. 30, 67-76, 1979) which has been extensively tested against experimental data and can be described by S = k • Cnt ~. The assumed sorption model allows parameters to be obtained by independent batch and column experiments. A numerical technique is used to solve the solute transport equation incorporating a correction to numerical dispersion to improve the numerical solution. An analytical solution for a simplified case is also presented to test the numerical technique. Parameter sensitivity analysis shows that influent concentration and the parameter k strongly affect the initial breakthrough time of solute, with m and n affecting the shape of the breakthrough curve. Preliminary investigations show that the applicability of the model to describe column experimental breakthrough curves is promising. Key words--phosphorus transport, time-dependent sorption, two-consecutive reactions, solute transport
modelling
INTRODUCTION
Understanding the movement of inorganic phosphates in soil is essential to predict the potential of phosphorus pollution of groundwater from fertilizer application. Management of the problem can be greatly assisted if a model is available that can quantitatively describe the movement of phosphate with infiltrating water and its interactions with soil. Phosphorus contamination of groundwater from application of inorganic and organic fertilizers has in many cases resulted in euthrophication of receiving waters when the groundwater discharges to shallow surface water bodies. The dangers of excessive phosphorus leaching are well documented (Birch, 1982; Lukatelich et al., 1987). Reactions of phosphates with soils are time-dependent (Barrow and Shaw, 1974; Munns and Fox, 1976; Ryden et al., 1977) and several models have been proposed. The reactions are generally viewed as either phosphate adsorption by or resorption from the soils. Selim et al. (1976a) combined nonlinear-kinetic and instantaneous terms, Cameron and Klute (1977) developed a model using a combination of linear-kinetic and linear-instantaneous terms. These models assume that soil particles consist of two different reaction sites. De Camargo et al. (1979) also used the two-site model. Their model was similar to Cameron and Klute's (1977) model, but they assumed that the kinetic-type reaction was due to adsorption sites which were less accessible. Mansell et al. (1985) modelled phosphate transport in a sandy soil by using
the three compartment model of Barrow and Shaw (1975). The model assumes that sorption-desoprtion follows two consecutive steps. During the first step, physical adsorption rapidly removes phosphate from solution, and during the second step the adsorbed phosphate undergoes chemisorption following an n th order kinetic reaction. Recently, Gerritse (1989) derived a phosphate transport model for sandy soils incorporating a time-dependent reaction. The reaction of phosphate and soil in his model was assumed to proceed according to the first-order rate equation of Munns and Fox (1976). The above models have been tested against experimental data and are potentially applicable to describe phosphate transport in soils. Generally, the models include kinetic terms with parameters which are difficult to obtain independently. In the models proposed by Selim et al. (1976a), Cameron and Klute (1977) and De Camargo et al. (1979), the proportion of the two different sites to produce a good fit of the model to the observed breakthrough experimental data has to be obtained by trial and error, while for the models of Mansel et al. (1985) and Gerritse (1989), the experimental breakthrough curve is needed to determine the ratio of forward and backward physical adsorption. It would be desirable to develop a solute transport model which allows an independent evaluation of the sorption-desorption parameters. Although there have been many attempts to model the time-dependent sorption of phosphate by soils, the model proposed by Barrow and Shaw (1979) appears to have a theoretical basis and fits
1205
1206
SuPra NOTODARMOJOet al.
experimental data of phosphate sorption by soils from many parts of the world (Probert, 1983; Barrow, 1989). No transport model appears to have been developed incorporating this more universally applicable sorption equation. The objective of this paper is to simulate the movement of phosphate in soils and groundwater in one-dimensional steady flow conditions, incorporating concentration and time dependent reactions, and with model parameters that can be evaluated independently from batch and column experiments. THEORETICAL
The two-consecutive reaction model o f P sorptiondes.orption
Sorption experiments are often conducted over a short period of time and equilibrium is assumed between sorbed and liquid phases. For phosphorus, true equilibrium is not reached within a short period. The reaction of phosphate with soils is complicated and the process is still puzzling. Much research has been undertaken to elucidate the processes involved. A plausible and likely mechanism has been suggested by Barrow and Shaw (1975) consisting of a threecompartment model for phosphate reaction with soils: A ~--,B ~--~C. Compartment A contains phosphate in the soil solution. Its concentration is determined by the concentration of adsorbed phosphate in compartment B according to an adsorption process. Compartment C contains phosphate which has been converted into a form which is not in direct equilibrium with compartment A. The reaction of dissolved phosphate (in compartment A) with the soil is considered as a rapid surface adsorption (or desorption) onto reactive surfaces followed by a slow penetration (or release) by solid state diffusion of this surface phosphate into subsurface horizons within the interior of the particles (compartment C). Barrow and Shaw (1979) concluded that the best description of the combined effect of both on phosphate sorption is of the form: S =k.C".t m
(1)
where S is the phosphate adsorbed onto the soil particles (~g/g of soil), C is solute concentration in the pore space (#g/cm3), k is a constant, which gives a relative measure of the adsorbent capacity of the soil for phosphate, n and m are constants which indicate the dependence of sorption on concentration and time (t) with a high value indicating stronger dependence. Barrow (1983) showed this simple descriptive model to be in good agreement with a more complex mechanistic model. A plot of log C vs log S from a series of batch experiments conducted at various concentrations and times of incubation would give the constants k, m and
n. The constants are affected by many factors such as: properties of the soils, pH, temperature, electrolyte background and soil-to-solution ratio (Barrow and Shaw, 1979; Syers and Iskandar, 1981). Phosphorus transport with two-consecutive reactions
In many cases the horizontal movement of pollutants such as phosphate in soils and groundwater can be simulated by a one-dimensional model (Selim et al., 1976b; Kinzelbach, 1986): [ O C ] ~ q C O S OOCst =--OzO Ov.O .~z - - ~ z - p "-~
(2)
where q is the volumetric flux of water through unit cross sectional area (cm/day), v = q/O is the average pore water velocity (cm/day) and Dv is the hydrodynamic dispersion coefficient (crn2/day). S is the solute associated with the solid phase (#g/g of solid), p is the bulk density of solid, C is the solute concentration in the solution phase (ktg/cm ~) and 0 is the volumetric water content. For steady uniform water movement the variables q, 0, v and D v are constant. Equation (2) then can be simplified to: OC 8t
OzC ~C D~. sz 2 - v c3z
p ~3S 0 Ot
(3)
Assuming S can be described by equation (1), then: 8S k . C , m.t,,_~ + k . t , , . n . C , _ l O C " dt = Ot
(4)
Substitution of equation (4) into equation (3) results in: OC c~2C - - =Dv.~z2 - V ~t
OC
"~ --~P . k . C " . m . t -o.[k'tm.n'C"-
m-l
i
8C
1~-.
(5)
Letting R = l + K ' t " ' n - C n ] and K = p / 0 " k , and rearranging equation (5) we then get OC
Dv 02C
O---t= R
dz z
v ~C
K
R Oz
R C"'m'tm-I"
(6)
Equation (6) is the solute flow equation with time-dependent sorption assuming no stagnant (immobile) water phase within the pores of the media. Analytical solution for a simplified case
Because of its non-linear nature, equation (6) is difficult to solve analytically. However, simplifying the equation allows analytical solution. When the value of parameters rn and n are equal to 1, equation (4) becomes: ~S c3t
~C 0t
-- =k.C +k.t.--
(7)
and equation (6) becomes: OC 02C 0C (1 + K . t ) - ~ = O~ . -~-fz2 -- v .-~z - K . C
(8)
Phosphorus transport in soils and groundwater Assuming a semi-infinite column length, a well mixed feed solution at the top boundary (z = O) and initially no sorbed or liquid phase phosphorus in the column, the appropriate boundary and initial conditions become: C=Co
at
z--O
~C dz
----0
as
C=S=0
at
t--0
for allz.
dC O2C OC K ' z "-~z = D "-~z 2 - V "-~z
(9)
-
K "C
(10)
with transformed initial condition: T=I.
(11)
Further using the Euler transform for z, z=C,
then
dz=e".dt?
(12)
and setting C(z, z) = ok(z, rl)
(13)
gives: K . ~ n~
C(z,t)
- 02~b 0~b = uv.-~z: -- V.~z - K . q b
at
t?=0;
a¢
--~0 Oz
q~=C0
as
f F R.z-v.t
z--,~.
at
(14)
z=0
Co f
(17)
Numerical solution for general cases When analytical solution is not possible, numerical approximation can often be used successfully. For one-dimensional solute transport, finite difference techniques are widely used. The explicit finite difference technique is commonly used, even though it may require extended computing time because of its restrictive stability criteria. However, by introducing corrections for numerical dispersion (Appendix A), the stability criteria become less restrictive. The (z, t) plane is discretized by a rectangular grid of points: i = 1, 2, 3 . . . . , I along the z-axis, and j = 0, 1, 2, 3 . . . along the t-axis.
(15)
The distance between the nodes Az remains constant:
Equation (14) together with the initial and boundary conditions [equation (15)] has been solved by van Genuehten (1981) to give the analytical solution: C(z,t)=T'~exP~
q
Figure 1 illustrates concentration breakthrough curves calculated by using equation (16) for various values of K. The number of pore volumes is defined as V/Vo where V is the total volume of liquid phase which has flowed through the column at a given time and V0 is the volume of liquid phase in the column. It can be observed that the curves exhibit a fairly flat top because of the strong dependence of sorption on time. A high value of parameter K gives a low final quasi-steady outlet concentration and a much delayed breakthrough.
with the following initial and boundary conditions: q~=0
i
Co = .er
2-7(B.-~Tt-)~,2 ].
z--.oo
at
column is available (Lapidus and Amundson, 1952):
for a l l t
Now if we let z = 1 + K . t , then d z / K = d t and we can rewrite (8) as:
C=0
1207
/(v
-
Az = L / ( I - 1)
~)z'x
2-~.-~ )
100--
x erfc; (K.z - ,,,,) \ 2 ' ~ 1
dT--
(18)
+ e x p/(v t--t \ 2-D, ,/
/j
x erfc( ( K . z + #~1) ~
tO
80-
0 tO X O3
(16)
K=0.2 K=0.4 K = 0.6K = 0.8
,f K=I K=2
D v : 25 cm2/d 1O0 cm/d
60--
t13 t~ p.
where #=v.
.o
1
+ 4 " K ' D , ~ 1/2 v----T--- ]
40--
cO3 °E -D
20--
rl = ln(l + K. t) and effc(x) is the complementary error function. With the sorption-desorption simplified to follow a linear isotherm, S = k . C, an analytical solution to the solute transport equation for a concentrationtype inlet boundary condition and semi-infinite
0.0
I 0.5
1.0
1.5
2.0
2.5
3.0
N u m b e r of pore volumes
Fig. I. The effect of parameter K on the breakthrough curves calculated with equation (16).
SUPRi NOTODARMOJOet al.
1208
where L is the column length and I is the total number of the nodes along the column. The time difference At is determined by the stability criteria of the scheme. The finite difference equation is derived by using a Taylor expansion (Mickley et al., 1957; Carnahan et al., 1969). Here we use a forward difference in time, a central difference for the second derivative, the dispersive term, and a backward difference in space for the advective term: 0C ~t
Ci +l - C ! At
Dv [C~+a - 2"C~ + CJ-I? R Az 2 e ~C~-C~:_,~ K t ~(19) R" L _ R" [Cj],.m . 1 where
R = l + K'tm'n'[C~]"
~.
(20)
Equations (19 and (20) together with appropriate boundary and initial conditions [equation (9)] can be used for the computation; however, a correction to the hydrodynamic dispersion is necessary to minimize the inherent differencing error and to improve the performance of the finite difference technique. Van Genuchten and Parker (1984) showed that with the concentration fixed at the inlet and a linear adsorption isotherm, a better result was gained if the
column was considered as a portion of an infinite column. The exit concentration is then simply the concentration at position z = L Figure 2 is a schematic diagram of the column system and the finite difference compartments.
Correction due to numerical dispersion Numerical approximation of a partial differential equation deviates from the exact solution. The deviation is a direct result of replacing the original differential equation [e.g. equation (6)] by a difference equation with a finite time step (At) and layers of finite thickness (Az). This deviation can be quite substantial and is called numerical dispersion (van Genuchten and Wierenga, 1974). It was first quantified by Chaudari (1971) as resulting from the neglect of higher order terms in At or Az in the Taylor series approximation, van Genuchten and Wierenga (1974) and De Smedt et al. (1981) incorporated more terms, and derived techniques to reduce the errors caused by numerical dispersion. They showed how these terms could be found in the case of classical hydrodynamic dispersion, in both the case of steady water flow and the case of transient unsaturated water flow with the presence of immobile water. Only the important terms of order At (the lower orders) were taken into account. By using a similar technique, a correction to the numerical dispersion is derived in our non-linear case (Appendix A) and the result can be written in the following form:
D* = Dv - ~ .(Az--AtR'V)
(21)
Co (influent concentration)
where D* is an effective hydrodynamic dispersion coefficient. When used in place of Dr, it should minimize numerical dispersion. 1 Top
of the c o l u m n
,=~---- 2 ,=-- 3
Az = Depth increment
Stability of finite difference scheme One method to test the stability of a finite difference technique uses a matrix method (Smith, 1978); the technique was used by Binh (1981). Using this method, the difference equation for the partial difference equation [equation (19)] is applied to all compartments in the system. Applying the Gerschgorin circle theorem (Appendix B) to this matrix, we find the stability criteria for the time-step is: At ~<
,,=,-- L-1
L
Bottom of the c o l u m n , Effluent concentrations L+I are computed at this point. L+2 Bottom of imaginary column, assumed at infinity, where concentratlon gradient is zero C e (effluent concentration)
Fig. 2. Representation of the column for numerical solution.
1
2.Dv v K R.Az 2 t - ~ . A z + ~ . R C n ' m ' t " - I
(22)
Another constraint that has to be satisfied is the mass balance of the system. The mass balance criteria can be expressed by the Courant and Neumann criteria (Kinzelbach, 1986). The Courant criterion requires that: Az At ~<--. v
(23)
In physical terms, the criterion states that no more
Phosphorus transport in soils and groundwater pollutant mass can leave a segment via convection during the time interval [t + At] than is inside at the beginning of the time interval. The Neumann criterion, which is analogous to the Courant criterion, ensures that the concentration gradients cannot be reversed by dispersion-diffusion fluxes alone. It has the form: ~LZ 2
At
(24)
~ < -
2-D*"
AND DISCUSSION
Validation of numerical solution The finite difference method is said to converge if the discretization error approaches zero as the grid spacings Az and At tend to zero (Carnahan et al., 1969). When the stability criteria are satisfied we have found that our numerical solutions are convergent. To assess the accuracy of the numerical solution, we compare numerical results to those from analytical solutions for the simple cases. This approach was used by van Genuchten and Wierenga (1976) in solving equations for the convective--dispersive trans-
100 - -
v ,. L ,n ,, m -
D
port of reactive solute with exchange between mobile and stagnant water. Figure 3 shows a comparison between the numerical approximation and the analytical solution as given by equation (16), while Fig. 4 shows a comparison with the analytical solution derived by Lapidus and Amundson (1952). It can be observed that the agreement between numerical approximation and analytical solution is very good in both cases.
Model sensitivity analysis
For a given calculation cycle, the value of At obtained with equations (22) and (24) will vary for every compartment. The minimum calculated At for all compartments is then used as the criterion for the next step. For most cases, equation (22) has been found to be adequate to ensure stability. However, for solute transport where the convective flux is dominant or where the solute exhibits a very strong reaction with the media, a combination of equations (22), (23) and (24) needs to be used. RESULTS
1209
100 c m / d 2 0 cm 1 1
To evaluate the sensitivity of the solution of the phosphorus transport model to different model parameters, breakthrough curves are calculated for a 20 cm long soil column leached with 3 pore volumes of phosphate solution followed by solute free water at a constant flow-rate. The results of this simulation are presented in Figs 5-9. The influent concentration has a profound effect on the effluent concentration (Fig. 5). The retardation or "hold up" time is large for low influent concentrations, and for high concentrations, the effect of sorption is reduced and the shape of the c u r v e converges to that of a non-reactive solute case, with an S-shaped curve. When the incoming solute concentration is low, the gradient for transfer into the grains becomes low, hence the process of transfer to the solid occurs only at a very low rate. This does not lower the concentration at the surface and consequently lessens the surface reaction rate. This mechanism explains why at lower input concentrations applied to the soil, the desorption curve shows more tailing. When K increases, the hold up time also increases significantly although not in a linear proportion, and
D v" 5 cm2/d K - 0.5
NA
80-0
C 0
.--_ x u:,
--
:
60
NB
--
D v" 50 cm2/d K-2.5
u) (I)
-~
NA
0
"~ C
D v - 100cm2/d
4 0 --
K ,, 5 . 0
O'-E - NB~
20-
0
1
I
I
I
I
I
2
3
4
5
6
Number
of pore volumes
Fig. 3. A comparison between analytical solution and numerical solutions. AN = analytical solution, NA = numerical solution with correction and NB = numerical solution without correction.
1210
SUPRI NOTODARMOJOet al. 100
cm
~~'~=-~2
v = 100 cm/d
--
L = 20 n =l m-0
NA
v
80--
O tO O
N
/d
K=0.5
X
60--
Dv=5Ocm2/d K=1.5
t-
.2 tO
40--
E
~5 20--
NB 0
1
2
3
4
N u m b e r of p o r e v o l u m e s
Fig. 4. A comparison between analytical for n = 1 and m = 0 and numerical solutions. LA = analytical solution, NA = numerical solution with correction and NB = numerical solution without correction. the peak concentration decreases concurrently (Fig. 6). This is a direct consequence of the parameter K being a relative measure of the adsorbent capacity for phosphate and proportional to the reactive surface of the soils where the instantaneous reaction between soil particles and solute occurs. The magnitude of parameter n for the soilphosphate interaction is expected to be less than unity. The parameter n has a significant effect on the starting time of the breakthrough curves for low n values (Fig. 7). However, the slope of the
breakthrough and the descending portion of the curves are most affected by parameter n. The value of n relates to the availability of surface reaction sites and governs the ability of the particles to adsorb solute at higher concentrations. The slope of the curve is proportional to the ratio of parameter n to m. When n/m is close to unity, the shape of the curve is steeper since the rate of solid diffusion, which is responsible for the time dependency, is equal to the surface reaction rate at a given concentration.
Dv=25cm2/d v = 50 cm/d K=1.2 n=0.7 m=
100
8O O tO O
c 0
60
--
40
--
E
~5 20 --
0
1
I
I
I
,
2
3
4
5
Number of pore volumes Fig. 5. The effect of influent concentration on the breakthrough curves.
6
Phosphorus transport in soils and groundwater
1211
K =0.2
100
K=0.4 K =0.6 A
8O
K =0.8
o~ v
K=I.0
O ¢O O x (D tO tO O
.£
K=2.0
60
K=4.0
40
tO ri f
E
~5 20
I 0
I 1
2
3
4
5
6
Number of pore volumes Fig. 6. The effect of parameter K on the breakthrough curves. An increase in m lowers the peak of the curves (Fig. 8). The parameter m affects the potential of phosphate storage in the interior of soil particles; a high value increases this potential and means that more can be absorbed during the sorption period and more can be leached from the soil particles during the leaching process. The effect of dispersion is illustrated in Fig. 9. Dispersion affects the starting point of the breakthrough, with a relatively little effect on the shape of the curves.
~
100 - -
A
80--
o~ v 0 r0 0 .~_ x q~
The general pattern of the relative exit concentration distribution is a non-symmetrical curve, with a sharp front in the ascending portion, followed immediately by a decrease in effluent concentration after displacement with solute-free solution, then a significant tailing takes place. The net effect of a change in one parameter is not independent of the values of other parameters. The extent of tailing has important implications for the phosphate decontamination of soils and groundwater. Once contaminated by excessive phosphates, a
60--
n=l.0 n=0.9 n=0.8 n=0.7 n=0.6 n = 0.5
to
C
40
C
E ~5
20
0
1
2
I 3
4
N u m b e r of p o r e v o l u m e s
Fig. 7. The effect of parameter n on the breakthrough curves.
5
6
1212
SUPRI
100 -
NOTOD.CtMOJOet al.
D v = 25 cm 2/d v = 50 cm/d K=1.2 C O = 10 mg/I n=0.7
80 A
o~ 0 0 t.)
60--
x o') ~D
._(2
40--
r,
E ~5
20--
1
2
3
4
5
6
N u m b e r of p o r e v o l u m e s
Fig. 8. The effect of parameter m on the breakthrough curves. soil may take considerable time to desorb or displace sorbed phosphates due to the time-dependent sorption characteristics of phosphate in soil. Even with relatively low concentrations, the sensitivity analysis shows that the displacement of phosphate from soils may occur over an extended period. An example o f model application To illustrate the practical application of the model, two breakthrough curves from two different soil types are presented. Batch experiments were carried out to obtain estimates of the parameters k, n and m, according to Barrow and Shaw's (1979) method. Leaching experiments were carried out by using 20 cm soil columns following the method of Mansell et al. (1985). Hydrodynamic dispersion coefficients v = 100 cm/d K-1.2 C o = 10 mg/I n=0.7 rn = 0.4
100 -
80-
were also determined by using a similar miscible displacement experiment technique. Details of experiments and more extensive results will be presented in a separate paper. Figure l0 shows the breakthrough curve from an Australian sandy soil (Gavin A horizon), for which the important properties were given by McArthur and Bettenay (1960). Figure 11 shows the breakthough curve from an Indonesian loamy clay soil (Buah Batu soil). As can be observed, the agreement between experimental results and simulated curves are good. Further development We have developed a one-dimensional model for the transport of soluble phosphate incorporating time dependent sorption using a relationship that has been
~
Dv = 5 cmS/d
V. oc..
'
O tO 0
\\~--- D,=IO
60-
cm2/d
X (/)
C 0
4O
• Dv - 50 cm2/d
C
E
a
• Dv = 60 cma/d
2O
Dv
70
[
I
3
4
Cm 2 /dl
m o
2
~~
I 5
N u m b e r of p o r e v o l u m e s
Fig. 9. The effect of hydrodynamic dispersion on the breakthrough curves.
6
Phosphorus transport in soils and groundwater
A
100
--
80
--
1213
[ m
O e-" 0 O
D v " 15.50 cm2/d v - 6 3 . 1 5 cm/d Co - 2 . 5 0 mg/I K - 0.46 n = 0.46 m ,= 0 . 2 3
n
60
--
X
O
._o
40
--
20
--
t,--
E
a
O
cP []
I 0
1
I
I
I
I
2
3
4
5
Number
of p o r e v o l u m e s
Fig. 10. Experiment and simulated breakthrough curves from leaching with 2.50 mg/! phosphate solution for Gavin A soil.
tested over a wide range of soils and conditions. The model allows the independent determination of the sorption parameters. This is an improvement over similar existing models. A one-dimensional model can be applied in a wide range of field situations under restricted conditions where one-dimensional flow can be approximated. 100
Real systems are, however, generally very complex. The presence of immobile water cannot be ignored and a number of models referred to in the Introduction account for the time-dependent behaviour of sorption simply by its presence. Time-dependent sorption of phosphorus in soils occurs, however, even in batch systems. The presence of immobile water
m
cb
/° J
B
80
-
l
D v - 118.60 cm2/d v - 6 7 . 0 2 cm/d C O = 4 0 0 rng/I K= 197 n = 0.26 m=0.2
[]
oo
A v O cO O
m
60
-
40
--
X
._m C
O O
r0
t,-
O [] 0
E
a
20--
I 0
I
I
I
I
I
6
12
18
24
30
36
Number of pore volumes Fig. 11. Experiment and simulated breakthrough curves from leaching with 400 mg/l phosphate solution
for Buahhatu B soil.
1214
SUPRI NOTODARMOJOet al.
would add to this in a flow system, although this effect would be small when pore velocity is low. One-dimensional flow in the vadoze zone is characterized by varying soil water content. We are currently developing our model to incorporate the unsteady state condition due to irrigation and subsequent moisture redistribution on phosphorus transport. In particular we are investigating whether the steady state model can approximate the behaviour of the more complicated unsteady state model. Results of this study will be published separately.
CONCLUSION A non-linear phosphate transport model with two-consecutive reactions has been presented. The model can use independent batch experiments to obtain the value of the parameters that characterize the phosphate reaction with soils. An explicit finite difference code has been developed incorporating a correction for numerical dispersion. A new analytical solution for simplified conditions was derived and evaluated to aid validation of the finite difference code. Numerical results showed little numerical dispersion and were in close agreement with the new analytical solution and the analytical solution of a simple linear adsorption case. Parameter sensitivity analysis indicates that the model is very sensitive to the influent concentration (Co), and the parameter k, affecting the starting time of the breakthrough curve. Parameters n and Dv affect the start of the breakthrough curves although the effect of Dv is not as significant. The model also shows that at low applied phosphate concentrations, breakthrough may still occur, and sorbed phosphate can be desorbed by leaching. Further testing of the model against experimental data is planned Acknowledgement--We would like to thank International
Development Program of Australian Universities and Colleges (IDP) for providing a postgraduate scholarship to one of us (S.N.)
REFERENCES
Barrow N. J. (1983) A mechanistic model for describing the sorption and desorption of phosphate by soil. J. Soil Sci. 33, 733-750. Barrow N. J. (1989) The reaction of plant nutrients and pollutants with soil. Aust. J. Soil Res. 27, 475-492. Barrow N. J. and Shaw T. C. (1975) The slow reactions between soil and anions: 2. Effect of time and temperature on the decrease in phosphate concentration in the soil solution. Soil Sci. 119, 167-177. Barrow N. J. and Shaw T. C. (1979) Effects of solution and vigour of shaking on the rate of phosphate adsorption by soil. J. Soil Sci. 30, 67-76. Bear J. (1979) Hydraulics of Groundwater, pp. 225-276. McGraw-Hill, New York. Binh A. (1981) Neutralization of acid wastewater by limestone-bearing sand. Ph.D. thesis, Murdoch University, Western Australia.
Birch P. B. (1982) Phosphorus export from coastal plain drainage to the Peel-Harvey Estuarine system, Western Australia. Aust. J. Mar. Freshwat. Res. 33, 23-32. Cameron D. R. and Klute A. (1977) Convective~lispersive solute transport with a combined equilibrium and kinetic adsorption model. Wat. Resour. Res. 13, 183-188. Carnahan B., Luther H. A., and Wilkes J. O. (1969) Applied Numerical Methods, pp. 429-530. Wiley, New York. Chaudari N. M. (1971) An improved numerical technique for solving multidimensional miscible displacement equation. Soc. Petrol. Engng J. 11, 277-284. De Camargo A. O., Biggar J. W. and Nielsen D. R. (1979) Transport of inorganic phosphorus in an Alfisol. Soil Sci. Am. J. 43, 884-890. De Smedt F., Wierenga P. J. and van der Beken A. (1981) Theoretical and Experimental Studies of Solute Movement through Porous Media with Mobile and Immobile Water,
pp. 207-210. Vrije Universiteit Brussel, Belgium. van Genuchten M. Th. (1981) Non-equilibrium transport parameters from miscible displacement experiments. Research Report 119, USDA-ARS, Riverside, Calif. van Genuchten M. Th. and Parker J. C. (1984) Boundary conditions for displacement experiments through short laboratory soil columns. Soil Soc. Am. J. 48, 703-708. van Genuchten M. Th. and Wierenga P. J. (1974) Simulation of one-dimensional solute transfer in porous media. New Mexico Agric. Exp. Sm Bull. 628. van Genuchten M. Th. and Wierenga P. J. (1976) Numerical solution of the convective-dispersion with intra-aggregate diffusion and non-linear adsorption. In System Simulation in Water Resources (Edited by Vansteenkiste G. C.), pp. 275-292. North-Holland, New York. Gerritse R. G. (1989) Simulation of phosphate leaching in sandy soils. Aust. J. Soil Res. 27, 55~6. Kinzelbach W. (1986) Groundwater Modelling: an Introduction with Sample Programs in Basic, pp. 211-265. Elsevier, Amsterdam. Kreyzig E. (1983) Advance Engineering Mathematics, pp. 822-824. Wiley, New York. Lapidus L. and Amundson D. R. (1952) Mathematics of adsorption in beds. VI. The effect of longitudinal diffusion in ion exchange and chromatographic columns. J. Phys. Chem. 56, 984--988. Lukatelich R. J., Schofield N. J. and McComb A. J. (1987) Nutrient loading and macrophyte growth in Wilson Inlet, a bar-built south-west Australian estuary. Estuar. Coast. ShelfSci. 24, 141-165. Mansell R. S., McKenna P. J., Flaig E. and Hall M. (1985) Phosphate movement in columns of sandy soil from a wastewater-irrigated site. Soil Sci. 140, 594i8. McArthur W. M. and Bettenay E. (1960) The development and distribution of the soils of the Swan Coastal Plain, Western Australia. CSIRO Aust., Soil Publ. No. 16. Mickley H. S., Sherwood T. K. and Reed C. E. (1957) Applied Mathematics in Chemical Engineering. McGrawHill, New York. Munns D. N. and Fox R. L. (1976) The slow reaction which continues after phosphate adsorption: kinetic and equilibrium in some tropical soils. Soil Sci. Soc. Am. J. 40, 46-51. Probert M. E. (1983) The sorption of phosphate by soils. In Soils: an Australian View Point, pp. 427-435. CSIROAcademic Press, Sydney. Ryden J. C., McLaughlin J. R. and Syers J. K. (1977) Time-dependent sorption of phosphate by soils and hydrous ferric oxides J. Soil Sci. 28, 585-595. Selim H. M., Davidson J. M. and Mansell R. S. (1976a) Evaluation of a two-site adsorption-desorption model for describing transport in soil. Summer Computer Simulation Conference, 12-14 July 1976, Washington, D.C. Selim H. M., Mansell R. S. and Elzeftawy A. (1976b) Distribution of 2,4-D and water in soil during infiltration and redistribution. Soil Sci. 121, 176-183.
Phosphorus transport in soils and groundwater
1215
Smith G. D. (1978) Numerical Solution of Partial Differential Equations: Finite Differences Methods, 2nd edition. Oxford University Press, Oxford. Syers J. K. and Iskandar I. K. (1981) Soil-phosphorus chemistry. In Modeling Wastewater Renovation (Edited by Iskandar I. K.), pp. 575-599. Wiley, New York.
Further simplification by neglecting third order terms and non-linear terms in ac/at results in:
APPENDIX A
A comparison of equation (A7) with equation (6) suggests that an approximation to the effect of numerical dispersion from the forward time difference can be obtained when the difference equation uses an apparent or corrected dispersion coefficient: At .v 2 D*=Dv-~ 2.RZ. (A8)
Correction for Numerical Dispersion Consider the following Taylor series expansion of C{ + 1:
-C{+OC.Atm Cij+l _ Ot
02Catz 03CAt 3 + ~ T - - f -+ Ot3 6
+ terms with higher orders of At.
(AI)
Substitution of equation (6) into equation (A1) and rearrangement of the terms results in:
At
At
-
Oz2
R Oz
AC At - -
K =
----
v OC . C n . m . t m - l
R 3z
_[~
Oz2
O C
R Cn'm "tm-I
A C
Oz
J J C~+, Ci_, CiJ --
Az
2.Az
At 2 02 rDv 02C
v OC
K C".m.t m-I]
R az
R
AC At (A2)
R
v.Orl-]Oc
ot.Oz2
ot K
v O2C ROtOz
,
O [ ~ ' C "m'tm-I].
C{+l_C{
At
az 2
(A9)
-
1 '
~ t --K'[C~]"'m't"-t
_ vIc{-~c{_,
The second term on the right-hand side of equation (A2) can be considered as a first estimate of the effect of numerical dispersion (Chaudari, 1971; van Genuchten and Wierenga, 1974; De Smedt et al., 1981). Note that this only considers the effect of the approximation to numerical dispersion in time. This term is:
Dv 03C
J Ci_,
Az
+2
+ terms with higher order.
OF1 ] o 2 c
--
Az -C~+, - 2C~ + C~_ 1]
R Oz
+-6-~?[~ 0 :
At v2]02C
+T~J~z~ (A7)
It is known that central differences best approximate the spatial derivatives (van Genuchten and Wierenga, 1974). However, with backward differences in the advective term, equation (AS) has to be modified. Here we approximate the derivative with central differences and look for the error
AtO [ ~ 02C vOC K . c . . m . t m _ , ] + 2fit
- - - - - -
R
x
[C~+ t - 2C~ +
(AI0)
Comparing equations (6) and (AI0) the numerical model should use the following corrected hydrodynamic dispersion coefficient: (A3)
dt
van Genuchten and Wierenga (1974) and De Smedt et al. (1981) proceed as follows. As a first order approximation, only the advective term in equation (A3) is retained.
APPENDIX B
Stability Criteria of the Finite Difference Equation _ _
~ _ _
Ot
C n . m .l m-I
=
At
Oz2
R Oz
R
At O F
+-~
v OC]
L-~zj.
(A4)
Rearranging terms and expanding the last term, gives:
AC At
-
-
=
K C,.m.tm_ t R
--__.
v dC R Oz
At v 0 0C 2 R Oz Ot
DvO2C R c3z2 At OC O F l-]
AC -
=
At
--__
K.c,.m.tm-I R
Uj = A • U j - ]
(A5)
Considering equation (6) for OClOt and neglecting the last term of equation (A5): -
Finite-difference approximations contain errors which occur due to neglecting terms of higher order. These can often give rise to numerical inaccuracy or regions of solution instability. The stability of a finite difference scheme can be assessed using a matrix method proposed by Smith (1978). The finite difference representation of the governing partial differential equation (using backward differences for the advective term), with the top and bottom boundary conditions, is given in matrix form as:
v OC DvO2C R Oz2
(B1)
where U s is the vector of all values of Us at time j, U j is the vector at time j - 1, and A is the matrix containing all constants. By recursion technique, equation (B1) can be expressed as:
R 0z
UJ = A .U j-~ = A .(.4 . U j-2)
At v D~33C At v2 O2C Atv2 0C O [ 1 ] 2 R R Oz3 ~ 2 R 20z 2 + 2 R OzOz
(A6)
....
A ' A . . I A ( A . U °) (B2)
where U ° is the vector of all values at time j = 0. Then equation (B2) becomes:
U s = (A) j. U°.
(B3)
1216
S U P R I NOTODARMOJO et
which has the general form:
Defining the error vector e = U - U*, where U* is the vector of exact values, then: ~. = U j -
U * = ( A ) J ' ( U O - U * ) = AJ'Eo
C j÷ ~= A • CL
I/Jl ~< 1, k = 1,2,3 . . . . . n
Application of Gerschgorin's circle theorem to matrix A gives:
lit - ( 1 - 2 D - V - K ) I < ( D ) + ( D
Az
j
(BS)
(I-2D-
V-K)-(2D+
(1 - 4 D - 2 V - - K ) ~ #
D = ~-~7k'
v-
v'At
az.R'
<~1 - K .
(BII)
For stability: - I ~ - 1 only when: (1 - 4 D - 2 V - K ) > ~ I
(B6) or
Defining the following terms:
Dr'At
V)
or
J I
(B10)
~
C{_,l
R . [ C ]i , - , C ~.m.t,~
+ V)
where (1 - 2 D - V - K) are the diagonal elements (A,.,), and (D) + (D + V) is the sum o f moduli (2,) of the element along the row excluding the main diagonal for segments other than the boundary. Therefore,
where/~ is the eigenvalue of A (n x n). Now the problem is confined to testing the eigenvalues of matrix A. One method of finding the bounds of eigenvalues is Greschgorin's circle theorem (Smith, 1978; Kreyzig, 1983). This theorem states that if 2, is the sum of the moduli of the element along the s row excluding the diagonal element As.s, then each eigenvalue of A lies inside or on the boundary of a least one of the circles IV - As.~l = 2~. The finite difference equation (19) can be rewritten as:
R [_
(B9)
(]34)
where E0 is the vector of round-off error from the beginning. According to Smith (1978), a finite difference scheme is stable only when ~, remains bounded as j increases indefinitely. For this to happen the eigenvalues of A must have modulus values less than or equal to unity:
c~C _ C,J+' - C { _ D v f C { + , - 2 " C ~ + Ot At R L Az2
al.
4D + 2V + K ~< 2.
K=K'At[c~],_].m.t,,_l
Substituting the values of D, K and K into equation (B12) results in:
R
equation (B6) becomes: At ~<
c ? ' = [D + v]. c~_~ +[1-2.D-V-K].C~+DCJ,+~.
(BI2)
1
2.D v
v
(B13)
K
RAz E"FR~-'~z +--2"R "Cn-I "m "t m-I
(B7)
In matrix form this is: CI
c~ c~
)+1
1
CI
(D+V)
(1-2D-
( D + V)
V-K)
(I-2D-
c:
D
V-K)
D
(B8) w
m
c~ . Cn+
( D + V)
(I-2D-
( D + V)
V-K)
D
(1-2D-V-K)
c~ _Cn+l