An alternative approach to calculating chemical equilibrium in soil reactions

Ecological Modelling, 20 (1983) 165 173 Elsevier Science Publishers B.V., Amsterdam - Printed in The Netherlands

165

AN ALTERNATIVE APPROACH TO CALCULATING CHEMICAL EQUILIBRIUM IN SOIL REACTIONS

ERNESTO BOSATTA

Swedish University of Agricultural Sciences, Section of Systems Ecology, Post Box, S-750 07, Uppsala (Sweden) (Accepted for publication 29 March 1983)

ABSTRACT Bosatta, E., 1983. An alternative approach to calculating chemical equilibrium in soil reactions. Ecol. Modelling, 20:165 173. In most simulation models dealing with complex chemical reactions in soils, the hypothesis of instantaneous equilibrium (IE) is adopted. The use of this hypothesis (a form of the adiabatic elimination hypothesis) is based partially on the fact that the relaxation times of these reactions are much shorter than the corresponding ones of the external forces, and partially on the current ignorance about the kinetics of most of these reactions, which results in assuming that the system moves forward in time through steps of chemical equilibrium. The alternative method proposed here, on the contrary, does not move stepwise through equilibrium values, thus avoiding the intermediate numerical complications. It makes use of a single scale factor which can be adjusted in order to get "as close as desired" to the concentration values generated with the conventional IE method. Furthermore, a practical criterium is also deduced in order to define a convenient order of magnitude for the scale factor.

INTRODUCTION S i m u l a t i o n m o d e l s dealing with a wide range of chemical reactions are n o w a d a y s b e c o m i n g m o r e frequent, since the negative effects of a m b i e n t a l p o l l u t i o n are generally m a n i f e s t e d t h r o u g h a long c h a i n of processes involving m a n y c h e m i c a l reactions. In a m o d e l a i m e d at c o n s i d e r i n g the effects of rainfall acidity in the leaching o f nutrients f r o m soils, Reuss (1980) c o n s i d e r e d i n t e r a c t i o n s bet w e e n several reactions: sulfate a d s o r p t i o n , c a t i o n e x c h a n g e involving Ca, Mg, H, A1, the d i s s o l u t i o n of gibsite, the several reactions involving C O 2 a n d c a r b o n a t e s in the soil, etc. Since the soil is an o p e n system, the d y n a m i c s of the r e a c t a n t s is n o t o n l y d e t e r m i n e d by the kinetics of the i n t e r n a l - - c h e m i 0304-3800/83/$03.00

© 1983 Elsevier Science Publishers B.V.

166 c a l - - r e a c t i o n s but also by external processes such as transport in the soil profile, weathering, etc. Information about parameters involved in external processes (such as the waterflow-rate) are, or can be made available, but in most cases even the basic kinetics of the chemical reactions are totally unknown. On the other hand, information concerning the parameters in the states of chemical equilibrium of the reactions is usually available as, for example, the equilibrium constants. The following procedure is the one commonly used in building models used to simulate the time dynamics of the concentrations of the reactants in the soil. The values of these concentrations are given at the beginning and are used to calculate the rates of external processes (leaching, for example). Next, these rates are used to update the values of certain variables--total variables--defined in such a way that their rates of change do not depend explicitly on the internal kinetics--an example of these variables would be the total amount of Ca, viz, adsorbed, plus in solution, if only exchange processes were considered. To move forward in time, a new value of the external rates is required, which in turn requires a new value for the concentrations. These values are calculated from the equations of chemical equilibrium (of the "law of mass action" kind) which are formulated in a way relating the concentrations to the total variables. This procedure is mainly based on the assumption that the internal processes are much faster than the external ones and that, in this way, chemical equilibrium is totally established at each time step. The main difficulty with this IE (instantaneous equilibrium) approach is the complexity involved in the last step, i.e., in solving the values of the concentrations fulfilling the equations of equilibrium, a procedure generally requiring the use of rather extended algor i t h m s - s u b r o u t i n e s - w h i c h are to be solved at each step of the integration program (see for example Reuss, 1980; Christophersen et al., 1982; Chen et al., 1982). As a consequence, the complexity of the model is greatly increased and an easy overview is lost. On the other hand, the question remains open whether so much numerical sophistication is justified when so little is known about the general virtues of the model. The aim of this paper is to present an alternative method of modelling the chemical equilibrium approach that avoids the use of intermediate algorithms. It is based on a formulation which mimics an approximate kinetics (AK) for the chemical reactions and makes use of a single parameter. Furthermore, the value of this parameter can normally be chosen in such a w a y that produced A K results lie as close as desired to the corresponding results of the IE method.

167 THE AK METHOD Let us consider a general system in which a number n of reactants are interacting through a number m (n > m) of chemical reactions and let C ~u) = (C~C2... Cn) denote the vector of concentrations in the soil of these reactants (these concentrations are not only limited to the bulk solution, they can also express concentrations in any arbitrary phase, the exchangeable phase being one such example). Let A be a particular reactant which takes part in several of the reactions. Using A as a reference we define the kinetics of these reactions in the following way. For each particular reaction, an equilibrium equation exists relating the concentrations of the reactants. Let C ~r) denote the vector of the actual concentrations of the reactants interacting with A in the c~th reaction. From the equilibrium equation we obtain CA = h r ( K r ,

C '~))

(1)

where K~ is the equilibrium constant and CAr is the concentration of A fulfilling the equilibrium condition for this reaction. Denoting with CA the actual concentration of A in solution, we define the "kinetics" of reaction c~ as proportional to the distance ¢,, = C A r - CA

(2)

In other words, the number of moles of reactant i transforming per unit time due to reaction ct is given by gia = k r P i r e a

i = 1 ... n a= l...m

(3)

where k~ (time ~) is a, for the moment, arbitrary specific rate, and v;~ is the stoichiometric number of i in reaction a. The total kinetics for reactant i is g,=Ek,d%%

(4)

where the summation extends over all reactions in which i participates. The dynamics of each reactant also depends on some external force, 4, (concentration, time-~), which is a function of c~N); SO, for reactant i

Using eqs. 4 and 5, the rate of change of C; is

But, as pointed out before, the specific rates, kr, are undefined. To make the model according to the A K method we assume first that k~ = k for all a,

168

so eq. 6 is transformed to

where k is now a common parameter for all the equations. In this way, a system of n differential equations is set up for the n reactants which can be solved in a straightforward way by any numerical method of integration (see for example, Lohammar, 1979) provided that the initial conditions (concentrations) are defined at time t = 0. A COMPARISON BETWEEN THE METHODS

We ask now if any relationship exists between the solutions of model 7 - - t h e A K m o d e l - - a n d the solution of the corresponding IE model and, in particular, in which way this relationship depends upon k. To answer this we first have to put eq. 7 into a different form. From the n-coordinates C, a set of n-m new coordinates X, are first generated by means of a suitable linear transformation. The new X are such that their time derivatives do not depend on any e, i.e., X has zero projection in the m-dimensional subspace defined by the m linearly independent vectors %. Keeping m of the old concentrations as the remaining coordinates, the A K model will look like

2 j = C ( x , c) C=+i(C)+k Z

c)

j = 1,2 ....

(8a)

,= 1,2 .... m

(8b)

~=1

where the Fj are now a linear combination of the old q~i, C is a vector of dimension m and ~,,~ are the elements of the stoichiometric m x m matrix r (notice that when defining ~i~ = 0 if reactant i does not take part in reaction e~, v becomes a square matrix). In the IE method one solves the equations

X j = F j ( X , C)

j = 1,2 . . . . n-m

(9a)

E , ( X , C) = 0

a = 1, 2 . . . . m

(9b)

where, as described before, Xj are the total variables updated at each time step and eq. 9b defines the m equilibrium conditions which must be fulfilled at each time step. This means that at each step of the integration, a system of m simultaneous algebraic equations in the m unknown C must be solved. Let us change once more the face of the A K equations, this time by multiplying eq. 8b by ~,- ~, i.e. the inverse of the stoichiometric matrix; in this way

169

Xj = Fj( X , C) ~',~=f~(C)+k~,~(X,

C)

j = 1, 2 .... n - m

(10a)

a = 1 , 2 .... m

(10b)

where Y = v ' C and, similarly, f = v-'q~. It can be seen now (to this end multiply eq. 10b by k - ' ) that by letting the scale factor k grow, eq. 10b will increasingly approach eq. 9b and, in particular, when k = oz the eqs. 10 of the A K model will coincide with the eqs. 9 of the IE model. Suppose that a steady state point (i.e., an isolated, asymptotically stable point) exists for the system of eqs. 9a, b. Let us represent this ss point of the IE model with the n-dimensional vector ( x ' E c I E ) . Evidently, this is not an ss point for the system 1 0 - - t h e A K m o d e l - - b u t , it can very nearly be if k is chosen "sufficiently large". We assume then that an ss point of the AK model, ( x A K c A K ) exists in the neighbourhood of ( x I E c I E ) and we are going to estimate the distance between them (for a more rigorous and general treatment of this sort of problem see Tikhonov et al., 1970; G 6 b b e r and Seeling, 1975). To this end, the eqs. 10a, b are linearized around the point (x~Ec~E):

k-~f',~--k-~f,,(cIE)+

Y'~ B,~id ~ a = 1 , 2 .... m i

(11)

1

where Oo~i=

~ii

Ct = ,~tII

(12)

and d i = C i -- Ci IE

(13)

is the "distance" vector. At the ss point (xAKcAK), Y = 0, and so, from eq. 11 d~-k

' B - l f ( C 'E)

(14)

where B ~ is the inverse of the matrix defined in eq. 12. The norm of this vector, d, can now be used as an estimation of the distance between the ss points. Equation 14 shows that this distance decreases by increasing k. But, the important, practical point about this equation is that it relates k with the external forces, thereby providing a criterion to choose the order of magnitude of k; namely, if % is the shortest relaxation time included in the external forces, k should be chosen in such a way that k >> ~-~1.

170 A SIMPLE ILLUSTRATION

The previous ideas are illustrated here through a very simple case. Let us have A, B with concentrations CA and C B and the chemical reaction

A~B

(15)

with the equilibrium law

K = CA/C a

(16)

Selecting A as a reference we define the kinetics of reaction as proportional to =

KC

B -

(17)

CA

Furthermore, the stoichiometric numbers are ~'a = 1 and u a = - 1. As external forces we have t w o - - c o n s t a n t - - r a t e s of inputs i A and i~ and two rates of leaching XCA and hC B (for simplicity it is assumed that A and B have the same specific rate of leaching, h). So, according to eq. 7, the A K model is given by

#,, = i , , - XC,, + k ,

(lSa)

C"a = i B - XC B - k~

(lgb)

where k is the scale factor. We now define the " t o t a l variable" X = CA + C B

(19)

and we keep CA as the remaining variable. In these new variables the A K model is ) ( = i A + i R - X X = i - XX

CA=iA-hCA+k[KX-(K+

(20a) I)CA]

(20b)

which are the equivalents of eqs. 8 or 10 (eqs. 8 and 10 are coincident for this problem since only one reaction is involved). If instead, one thinks in the IE mode, the following eqs. would be set up for this problem:

X=i-XX

(21a)

K X - ( K + 1)CA = 0

(21b)

which are the equivalent of eqs. 9a, b. Since this particular problem can be solved exactly, the whole trajectories generated with both approaches can be compared. It is verified that the distance between them decreases with increasing k. Here, we shall concentrate around the steady state points.

171

For the IE model, this is found from eqs. 21 as X

TM

(22a)

= i/~

CAIE=

K K+I

i X

(22b)

The ss point for the A K model is found from eqs. 20 (23a)

X AK = i / X cA K =

(23b)

iA + kKi/X

X+k(K+

1)

Now, according to a previous discussion (eq. 14), we choose k in such a way that k >> X. Introducing this approximation we get, using eqs. 22b and 23b i c 2 K - c']EI-=

iK

1):

(24)

as an estimation of the distance between the ss points. This value can also be derived directly from eq. 14 where, for this particular problem, the matrix B-1 reduces (see eq. 20b) to B-l=

(K+

1) -1

(25)

Equation 24 shows that in the limit k = oc both methods will provide the same steady state values. DISCUSSION

The hypothesis of instantaneous equilibrium (IE) is nothing but a form of the adiabatic elimination hypothesis, which can be applied to find approximate solutions to problems in which a set of parameters changes "slowly" in time compared to the changes of another set of variables (Haken, 1977). In this way, elimination of fast relaxing variables from a problem (reducing its complexity), can be a very useful tool in helping to analyse it (Bosatta, 1981; Bosatta and Staaf, 1982). What is questioned here is the usefulness of the same approximation when used in its purely numerical form. In its simpler form, and if m chemical reactions are involved, the IE method will require the inversion of a square matrix (the buffer matrix) at each step of the integration program (this matrix is the one defined in eq. 12 and is derived for the IE eqs. 9b if this algebraic system is substituted by a linear system of equations (Bosatta, 1980)). The use of intermediate iteration programs of this kind increases the size (and the clumsiness) of the simulation model, making communication more difficult. On the other hand, the A K model can be set up in a straightforward way, requiring the same basic information as the IE, viz, the

172

m equilibrium equations (note that the A K model is composed of equations of the 7 or 18a, b types; all the remaining, rather complex transformations used here, are only necessary for comparison with the IE model). Actually, what I am looking for is a reproduction of the results of the IE method in a simpler way. In order to strictly get the same results one should take k = ~ (eq. 14). The problem with the A K method is that if the scale factor chosen, k, is too large, the time step in the integration program must be chosen too short, thereby increasing the time of run of the model. On the other hand, I have shown that by choosing k some order of magnitude larger than the specific rates of the external forces, a good approximation to the IE results can be attained. In problems concerning reactions in soils, the main external driving force is generally produced by the flow of water through the system and so the main relaxation time is associated with the rate of percolating water. Of course, during the progress of a simulation with the A K model it is always possible to check how far the simulated concentrations are from the values fulfilling the various chemical equilibria. Consequently, in a model dealing with several reactions (exchange, A1 dissolution, etc.) I found that it was enough to take k --- 5k (where k is a specific rate of leaching estimated from typical values of rates of percolation), in order to have discrepancies less than a few percent between the simulated concetrations and the corresponding equilibrium values. Is this necessary to a greater extent in this kind of model? In conclusion I have shown that, since the relaxation times of the hydrological processes are normally of a moderate size, differential equations such as eq. 7 or eq. 18 can be of practical use to simulate chemical processes in soils. ACKNOWLEDGMENTS

I wish to thank F. Andersson and G. ,~gren for reading and commenting on the manuscript. This study was started within the Swedish Coniferous Forest Project and funded by the Swedish Natural Science Research Council. REFERENCES Bosatta, E., 1980. Modelling of soil processes--an introduction. In: T. Persson (Editor), Structure and Function of Northern Coniferous Forest--An Ecosystem Study. Ecol. Bull. Stockholm, 32: 553-564. Bosatta, E., 1981. A qualitative analysis of the stability of the root-microorganism soil system. I. Carbon-nitrogen status and nitrogen mineralization. Ecol. Modelling, 13: 223-236.

173 Bosatta, E. and Staaf, H., 1982. The control of nitrogen turnover in forest litter. Oikos~ 39: 143-151. Chen, C.W., Gherini, S.A., Dean, J.D., Hudson, R.J.M. and Goldstein, R.A., 1982. Development and calibration of the integrated lake-watershed acidification study (ILWAS) model. Presented at ACS National Meeting, Las Vegas, NV, March 1982. Christophersen, N., Seip, H.M. and Wright, R.F., 1982. A model for streamwater chemistry at Birkenes, Norway. Water Resour. Res., 18: 977-996. G/Sbber, F. and Seeling, F.F., 1975. Conditions for the application of the steady-state approximation to systems of differential equations. J. Math. Biol., 2 : 7 9 86. Haken, H., 1977. Synergetics. An Introduction. Springer-Verlag, Berlin, 320 pp. Lohammar, T., 1979. SIMP--interactive mini-computer package for simulating dynamic and static models. In: S. Halldin (Editor), Comparison of Forest Water and Energy Exchange Models. ISEM, Copenhagen, 258 pp. Reuss, J.O., 1980. Simulation of soil nutrients losses resulting from rainfall acidity. Ecol. Modelling, 11: 15-38. Tikhonov, A.N., Vasileva, A.B. and Volosov, V.M., 1970. Ordinary differential equations. In: E. Roubine (Editor), Mathematics Applied to Physics. Springer-Verlag, Berlin, pp. 211.