UNSATCHEM: Unsaturated Water and Solute Transport Modelwith Equilibrium and Kinetic Chemistry
D. L. Suarez* and J. Simunek
ABSTRACT
Numerous models have been developed for predicting major ion
chemistry in the soil zone and in recharge to groundwater. Soils that
contain CaCO, are prevalent in arid and semiarid regions, as well as
in humid and temperate regions that have been glaciated or containcarbonate bedrock. Under these conditions, carbonate-solution reactions and ion exchange are the dominant chemical processes. In thismodel we couple one-dimensional unsaturated water and solute transport with a major ion chemistry routine and plant water uptake.
The model has several unique features, including expressions relating
reductions in hydraulic conductivity to chemical factors, predictionof CO; partial pressure in the root zone based on a CO, productionmultiphase transport submodel, kinetic expressions for silicate weathering, calcite precipitation-dissolution, and dolomite dissolution, representation of B adsorption using the constant capacitance model, a
new method for predicting cation-exchange selectivity, the option to
use Fitzer ion interaction expressions for high ionic strength, and aplant growth submodel that includes water, salinity, and O2 stress.
The chemical submodel considers equilibrium ion exchange, as wellas various equilibrium and kinetic expressions for precipitation anddissolution of soil minerals, including gypsum, Mg carbonates, andsepiolite. The use of a predictive submodel for CO, production andtransport allows for the calculation of CO, concentrations with depth
and time. This enables us to avoid the assumption of constant CO,
distribution or constant pH required by previous models. Use of
kinetic expressions for carbonate chemistry allows a more realistic
simulation of soil and groundwater solution composition as well as
simulations of carbonate redistribution and climatic change with time.

M

ODELING THE TRANSPORT and concentration of major soluble ions in and below the root zone is arequisite for predicting groundwater quality as well as
for managing irrigation and fertilization practices. Manymodels have been developed in the past 20 yr to quantify
the physical and chemical processes that affect the transport of major ions. The hydrological models, solutetransport models, and aqueous chemical equilibrium
models were developed independently and only later
were these three kinds of models coupled. Among these,
Jury et al. (1978) developed a model for steady-statewater flow in the unsaturated zone that considered ion
USDA-ARS, U.S. Salinity Lab., 450 W. Big Springs Rd., Riverside, CA 92507-4617. Received 2 Aug. 1996. *Corresponding author(dsuarez@ussl.ars.usda.gov)

Published in Soil Sci. Soc. Am. J. 61:1633-1646 (1997).

exchange and calcite equilibrium. Schulz and Reardon(1983) presented a mixing cell analytical model forgroundwater transport that considered steady-state water flow, and a simplified chemical model (no ion pairing
or complexation) with ion exchange and calcite equilibria. Forster and Gerke (1988) proposed a multicomponent transport model based on steady-state water fluxthat considered a mobile-immobile water concept andcarbonate equilibrium. Numerous other models havebeen developed for steady-state groundwater flow coupled with chemical equilibrium (Jennings et al., 1982;
Walsh et al., 1984; Miller and Benson, 1983; Narasimhan
et al., 1986; among others).
Simulation of soil water flow and soil chemical processes requires consideration of water and solute transport under variable water content conditions. Modelingof unsaturated water flow coupled to equilibrium chemistry has been undertaken by only a few researchers.
Robbins et al. (1980a,b) combined an unsaturated waterflow model with an ion exchange and carbonate chemistry and gypsum submodel that assumed fixed soil pH.
Subsequently, Wagenet and Hudson (1987) developed
a similar model that input fixed CO2 at various depths
instead of assuming fixed pH. Russo (1986) combined
the solution chemistry of the Robbins et al. (1980a)
model with the Bresler (1973) transport model. Yeh andTripathi (1989, 1991) coupled a generalized chemical
equilibrium model with unsaturated water flow. Application of the Yeh and Tripathi model to soil environments requires specification of total inorganic C andfixed pH. The model does not contain any provisions
for plant water extraction. Simunek and Suarez (1994)
developed a two-dimensional model with unsaturated
water flow and major ion chemistry.
Evaluation of these models under field conditions isvery limited. Only Robbins et al. (1980b) tested theirmodel by comparing its results with experimental dataobtained from a lysimeter study. However, in a fieldstudy, Dudley et al. (1981) indicated that the Robbins
et al. (1980b) model could predict salinity, but not specific ion concentrations. These difficulties may be related to heterogeneous water flow, as mentioned by the
researchers, but may also result from the assumptionsof chemical equilibrium and fixed pH and the mannerin which the calcite equilibrium is solved.
Most previous models assume chemical equilibrium

1634

SOIL SCI. SOC. AM. J., VOL. 61, NOVEMBER-DECEMBER 1997

between the solution and the solid phases and require
input of either CO2 or pH with depth, without provision
for time dependence or interaction with soil processes.
The requirement that CO2 be input leads to serious
limitations to the utility of the models, especially innear-surface environments where climatic factors suchas rainfall and temperature have an important influenceon soil CO2 concentrations. Several regression relations
have been developed to predict average CO2 from temperature, rainfall, or evapotranspiration. In addition tothese factors, there are other factors that affect soil CO2
concentrations, such^as soil hydraulic properties andporosity (Suarez and Simunek, 1993). Also, average values do not consider the important effect of seasonal andshort-term cycles at each site. Buyanovsky and Wagner (1983) measured seasonal changes in soil CO2 concentration of 0.3 to 8% for cropped fields in Missouri.
The CO2 model developed by Simunek and Suarez(1993) was able to predict CO2 production, transport,
and distribution in the soil based on water inputs, soil
hydraulic properties, potential evapotranspiration, anda production submodel.
The generally utilized assumption of mineral equilibrium is also only a rough approximation for soil andshallow groundwater environments. Suarez (1977a)
measured calcite supersaturation of groundwaters beneath two irrigated regions. In both regions, waters wereon average threefold supersaturated. Suarez et al. (1992)
also determined calcite supersaturation in soil water
sampled from the (water) unsaturated zone. The utilityof a kinetic approach to improve prediction of major ion
chemistry was shown by Suarez (1985), who presented
results of simulations using steady-state unsaturated water flow combined with calcite precipitation kinetics using the Plummer et al. (1978) rate equations.
Existing models do not consider the effects of chemical properties on hydraulic conductivity. It is well documented, however, that low electrolyte concentration,high sodicity (McNeal, 1968; Frenkel et al., 1978), andelevated pH (Suarez et al., 1984) all adversely affecthydraulic conductivity. Consideration of these factors
is essential to predicting water and solute movement insodic soils.
Our objective was to develop a one-dimensional unsaturated water flow and solute transport model for
predicting major ion and B chemistry in field environments, with emphasis on arid-zone chemistry, including
the processes of plant water uptake, root and plant
growth, and without the need for specifying the CO2or pH.

z is the spatial coordinate [L] (vertically upward), and Q is
the sink-source term [T"1], which represents the water uptake
rate by plant roots.
The unsaturated soil hydraulic properties are based on aset of closed-form equations (van Genuchten, 1980), using
the capillary model of Mualem (1976) and modified here toaccount for the effects of chemistry on these properties. The
soil water retention and hydraulic conductivity functions are
given byBS -e r6(/z) = 6r\ah\n)"'[2]andK(h) = rKsKr = r.

1 - (1 - Sl/m )"f

[3]

n >1

[4]

respectively, wherem = 1 - I/n„e - ere, - er[5]

and where 6r and 9S denote residual and saturated water content [L3 L~3], respectively, Ks is the saturated conductivity[L T"1], Kt is the relative hydraulic conductivity, r is the chemical reduction function, Sc is relative saturation, and m, n, anda. [LT1] are the empirical parameters of the hydraulic characteristics. Hydraulic characteristics are determined by a set ofsix parameters: 6r, 6S, a, n, r, and Ks.
The chemical composition of the soil water and the exchangeable Na percentage have a dramatic effect on the soil
hydraulic properties, as mentioned above. The hydraulic properties described by Ks and K, are optimal values, inasmuch as
the measurements are made under favorable chemical conditions. Numerous studies have documented that elevated levelsof exchangeable Na, especially in combination with low electrolyte concentrations, will cause clay dispersion and/or swelling and subsequent reduction in soil hydraulic conductivity.
In addition, Suarez et al. (1984) determined that elevatedlevels of pH have a very adverse effect on saturated hydraulicconductivity. The pH effects were determined from experiments in which the pH was varied at numerous fixed exchangeable Na and salinity levels. This allows separation of the pH
effects from the salinity and Na effects. The overall chemical
hydraulic reduction function, r, is given by

r =W

[6]

where r{ is the reduction due to salinity-sodicity interactions,and r2 is the reduction due to adverse pH. The r, values arecalculated for soils of California based on the experimental

work of McNeal (1968), where
= 1-

ve

[7]

where v and w are empirical factors, and e is defined by

MODEL DEVELOPMENT

6 = /Smec,i«e 3-6 X ICT2 INad*

Water Flow

where /smec,i,e is the weight fraction of smectite in the soil, d*
is an adjusted interlayer spacing, and *Na is the Na charge
fraction of the cation exchange. The reduction factor r2, representing the effect of pH on hydraulic conductivity, is calculated
from the experimental data of Suarez et al. (1984):

The Richards equation is used to describe one-dimensional
water movement in partially saturated incompressible porous
media:dt

dzt \dz[1]

where h is the water pressure head [L], 0 is the water content[L3 L~3], K is the hydraulic conductivity [L T"1], t is time [T],r2 = 1
r2 = 3.46 - 0.36pH
r2 = 0.1[8]

for pH < 6.83

for pH = 6.83 to 9.3
for pH > 9.3[9]

SUAREZ & SIMUNEK: UNSATCHEM MODEL FOR UNSATURATED TRANSPORT

In view of the differences among soils, these specific values
may not be generalized predictors of soil hydraulic conductivity but they do serve to illustrate the changes that occur ininfiltration and solute movement under various chemical conditions. At the present time we assume that the hydraulicchanges are reversible and that the water retention function
is not affected by the solution composition. Information onthese effects is not available at present.
Solute Transport Equation for Unsaturated

Porous MediaThe partial differential equation governing one-dimensional advective-dispersive chemical transport under transient
flow in partially saturated porous media is taken as
d6cT
dt

1635

The model considers both gas and aqueous-phase CO2transport and air-soil gas exchange and separate productionresponse functions for root respiration and microbial activity.
These response functions consider temperature, water, andO2 stress as well as plant growth. Heat transport and soil
temperature are calculated as described in Simunek and Suarez (1993) using air temperature as an input. The model canbe used for near-neutral to high pH environments as well as
for representation of the addition of acid to calcareous soils(alkalinity can be negative, where the acid will dissolve calcite
and pH increases quickly back to values where Al is notimportant). The model presently does not consider Al chemistry, thus is not suitable for simulations in which these reactions
are important.
Precipitation of Calcitede,.

The precipitation or dissolution of calcite in the presenceof CO2 can be described by the reaction

dt

i = 1, ns

dz

[10]

CaCO3 + CO2(g) + H2O ^ Ca2+ + 2HCOr[14]

where the solubility product K$p is given by

where CT, is the total dissolved concentration of the aqueouscomponent i [M L~3], CT, is the total sorbed concentration of
the aqueous component / [M M"1], CT, is the nonadsorbed
solid-phase concentration of aqueous component / [M M~'],
p is the bulk density of the soil [M L~3], D is the dispersion
coefficient [L2 T"1], q is the volumetric flux [L T"1], and «s is
the number of aqueous components. The coefficient D is thesum of the diffusion and dispersion components[11]
where T is the tortuosity factor, Dm is the coefficient of molecular diffusion [L2 T"1], and X is the dispersivity [L]. The water
flow, solute and CO2 transport equations are solved using theGalerkin finite element method, as described in Simunek andSuarez (1993).

K§r = (Ca2+)(COn[15]where parentheses denote activities. Substituting the equation
for the dissociation of water, Henry's law equation for solubility of CO2 in water, and the equations for the dissociation of
H2CO3 in water (Eq. [A2], [A3], and [A4] in the appendix)
we obtain:(Ca2+)(HCO3-)2 =^/7co2 (H20) =[16]

Rearranging and expressing in concentration yields the thirdorder equation[Ca2+ + x][HCO3- + 2;t]2 =[17]

Each of the nine major aqueous components (Ca, Mg, Na,
K, alkalinity, SO4, Cl, B, and Si) is defined by summation of
the individual species containing that component. For example
Ca and B are defined bywhere x is the quantity of calcite that must be dissolved or
added to the solution in order to attain equilibrium. Solutionof a second-order equation based on Eq. [15] will not converge
properly, requiring up to 1000 iterations, as CO2" is almost
always a very minor species in solution. This problem is present
in several currently utilized models.

CaT = [Ca2+] + [CaSOS] + [CaCO§]

Precipitation of Gypsum

Equilibrium Chemistry

BT = [B(OH)3] + [B(OH)4-][CaHCO3+][12]

where brackets represent concentrations. Alkalinity is definedbyAlkalinity = [HCO3-] + 2[COl~] + 2[CaCO3)]+ [CaHCO3+] + 2[MgCO§] + [MgHCO3+]+ 2[NaCO3~] + [NaHCO§] + [B(OH)4~]- [H+] + [OH-][13]
The reactions in the CO2-H2O system and complexation
reactions for major ions have been described in numerous
publications, thus further discussion is not needed. A list of
the variables and solids considered are given in Appendix A.
The CO2 partial pressure (pCO2) is input to the chemicalroutine after calculation by the CO2 production transport submodel. This subroutine is taken from the CO2 model describedin Simunek and Suarez (1993) and uses the parameters developed in Suarez and Simunek (1993).

The gypsum saturation status of a solution can be evaluatedby comparing the ion activity product, IAF, to the Ksf where
the IAF is defined by[Ca2+][SO|-] =

IAF

[18]

To obtain equilibrium, i.e., when the IAF is equal to thesolubility product K§r, a quantity of gypsum must be addedor removed from the solution. This quantity is obtained bysolving the quadratic equation.
Magnesium Precipitation
We consider that Mg precipitation can occur as a carbonate(either nesquehonite or hydromagnesite) or as a silicate (sepiolite). Although there is still controversy as to which of theMg carbonate phases is thermodynamically most stable, thisis not of concern in this evaluation. Since we are developinga predictive model, we consider only phases that either precipi-

1636

SOIL SCI. SOC. AM. J., VOL. 61, NOVEMBER-DECEMBER 1997

tate under earth surface conditions or occur frequently andare reactive under these conditions. Magnesite can thus beneglected, as it apparently does not form under earth surface
temperatures, is relatively rare, and its dissolution rate is exceedingly small. Similarly, we do not consider precipitationof dolomite, as true dolomite appears to very rarely form insoil environments. Also, if dolomite is present in the soil, we
do not assume that the soil solution will reach equilibriumwith that solid. The dissolution rate of dolomite is very slow
(Busenberg and Plummer, 1982), thus we incorporate kinetic
expressions for dissolution.
Removal of substantial amounts of Mg by calcite precipitation (high Mg calcite) is possible at high Mg/Ca ratios but is
not presently considered, as the precipitation kinetics, mineralstability, dependence of Mg substitution as a function of precipitation rate, and biological activity (which controls the Mg
substitution) are unknown. If nesquehonite or hydromagnesite
saturation is reached, the model will precipitate the predicted
Mg carbonate. The Mg carbonate precipitated, combined withcalcite precipitation, will probably represent the mixed Ca-Mg
precipitate called poorly ordered dolomite or protodolomite.
However, the resulting solution composition is much different
than that produced by simply forcing equilibrium with respect
to dolomite, as the model forms this mixed precipitate(calcite + Mg carbonate) under conditions of approximately
three orders of magnitude supersaturation with respect todolomite. This result is consistent with the high levels of dolomite supersaturation maintained in high Mg waters (Suarez,
1977, unpublished data). Precipitation (or dissolution, if present in the soil) of sepiolite is also considered by the model.
Sepiolite will readily precipitate into a solid with a KSP greaterthan that of well-crystallized sepiolite. Formation of this mineral requires high pH, high Mg concentrations, and low CO2
partial pressure.
At 25°C and at CO2 partial pressures above 10~327 kPa,
nesquehonite is stable relative to hydromagnesite. The precipitation or dissolution of nesquehonite (MgCO3-3H2O) in thepresence of CO2 is described by an equation similar to Eq.
[14] with the solubility product K™P defined by= (Mg2+)(COi-)(H20):[19]

Substituting Eq. [A2], [A3], and [A4] from the appendix, we
obtainPO,

We solve for equilibrium in a manner similar to that used
for calcite, with a third-order equation. The precipitation or
dissolution of hydromagnesite in the presence of CO2 is described by

Mg2Si307.5(OH)-3H20 + 4.5H20 + 4CO2(g) ^2Mg2+ + 3H4SiOS + 4HCO3-[24]

with the solubility product A|P defined by„ _ (Mg2+)2(H4SiQ^)3(OH-)4[25]

Substituting Eq. [Al], [A2], and [A3] from the appendix, we
obtain[26]
Relatively little information exists on the Si concentrations
in soil waters, especially in arid zones. In soil systems, Siconcentrations are controlled by dissolution (and possibly precipitation) of aluminosilicates and Si adsorption onto oxides
and aluminosilicates. As a result of these reactions, Si concentrations in soil solution follow a U-shaped curve with pH. The
data shown in Fig. 1 were obtained by reacting eight aridzone soils from the western USA with deionized water at asoil/water ratio of 1:3 (by weight). For each soil, a series ofsamples were adjusted with either 0.1 M NaOH or 0.1 M HQ
to achieve pH values in the range 7.5 to 11. Additional acidor base additions were made periodically to maintain pH.
After 14 d, final pH values were measured and samples were
filtered through O.l-jjun filters and Si analysis performed using
the heteropoly blue method (Clesceri et al., 1989). The secondorder regression equation

2SiO2 = 0.001(6.34 - 1.43pH + 0.0819pH2) [27]where SiO2 is the sum of all silica species expressed in moles
per liter, had an r value of 0.84 (Suarez, 1977b). This relationship probably provides only a rough estimate of Si concentrations, but we consider it acceptable because it is used only torestrain Mg concentrations at high levels of evapotranspiration, when Mg concentrations become very high at low CO2
and high pH.
In this model, we utilize Eq. [27] to predict the Si concentrations in solution for arid-land soils when we utilize the sepiolite
precipitation option. Freshly precipitated sepiolite has been
prepared in the laboratory at IAF values of 10~35 comparable
to the KW listed by Truesdell and Jones (1974), thus a kinetic
expression for precipitation is not required.
An additional option (for nonalkaline soils) is to consider
the Si concentration to be controlled by inputs from mineralweathering and concentrated only by processes of evapotranspiration. In this case we utilize kinetic expressions for theweathering of selected silicate minerals as discussed below.

Mg5(C03)4(OH)2-4H20 + 6C02(g) ^5Mg2+ + 10HCO3-[21]

with the solubility product K& defined by
KSHP = (Mg2+)5(COn4(OH-)2(H20)4[22]
Substituting Eq. [Al], [A2], [A3], and [A4] from the appendix,
we obtain(Mg2+)'(HC03-)'° =[23]

Again we solve for the amount to be precipitated or dissolvedin the manner described for calcite and nesquehonite.
The precipitation or dissolution of sepiolite in the presenceof CO2 is described by

Cation Exchange
Cation exchange between the solid phase and the solution
can be described by the Gapon equation (White and Zelazny,1986):
cf (cf[28]

where z, and Zj are the respective valences of species i and j
and the overscored concentrations are those of the exchanger
phase (in mojc kg"1 soil). We apply this equation for Ca,
Mg, Na, and K exchange, and a summation equation for theexchangeable ions:

CT = Ca2

Mg2

Na

KH

[29]

1637

SUAREZ & SIMUNEK: UNSATCHEM MODEL FOR UNSATURATED TRANSPORT

0.8

Table 1. Adjustable parameters for calculation of activity coefficients.
a'Species
6t
at

0.6

Mg+2+

———————I————————I————————1———————I——oo~
pH-Si0 2 Relationship for 8 Western Soils
°

5.5
4.0
3.5
5.0
5.4
5.4
6.0
4.0
9.0
3.5
5.4
4.0
2.5

Na

K*
S

SO5r

0.4

COi"

W

HCOf
CaHCO?+
MgHC0
3
+

0.2

8

9

H
OH-

II

10
pH

Fig. 1. Dissolved Si concentration as a function of pH for eight western
U.S. arid-zone soils. Soils were reacted at 1:3 soil/water ratios for
14 d.

0.165
0.20

5.0

Ca"

HjSiO,"H2SiO5~
B(OH)4"
CaCO§
CaSOS

0.075
0.015
-0.04
0
0
0
0
0
0
0

-0.5

-0.45*-0.63$

MgCO°,
MgSOS

as done previously in Robbins et al. (1980b). It is assumed
that the cation-exchange capacity CT is constant.
Existing models either require input of a set of soil-specificselectivity values or use a generalized value. We observe thatthe experimentally determined selectivity values are not constant across the range of exchanger composition for any given
soil and, further, that the variation from one soil to anothercan be large.
Recently, Suarez and Wood (1993) and Suarez et al. (1994)
developed predictive models for Ca-Mg and Ca-Na exchange,
respectively, that adequately predicted the ion selectivity fora range of smetitic soils. In their model, the selectivity value
is calculated from knowledge of the cation-exchange capacity(CEC) and organic matter content, and using the known mineral and organic matter selectivities, constructing a mixing
model for the overall selectivity. Initially the measured CEC
is partitioned into organic and inorganic components. If theorganic matter content is known, it is assumed that the CECof the organic matter (molc kg"1) is equal to its content in thesoil (kg kg"1) multiplied by its CEC. The calculations use avalue of 4.9 molc kg"1, the assumed CEC of pure organic
matter (Pratt, 1957). The clay CEC is obtained by subtracting
the organic matter CEC from the total. The resulting binary
selectivity values are represented by expressions such as(Mg2+)
L-^Mgclay + -''•MgomJ V^af~*T* /~~*

—

1 +

^Kom

[31]

1 + B VT

[32]

[33]

CECom =
CECc|ay-~*Kclay[34]

CECclay -1
1 + 1.96(Ca2+)
05/(Na+)J

CECc]ay = Jcaday + ^Naclay + ^Mgc|ay + ^Kday

[35]

[36]

[37]

[38]

where a' is an empirical parameter whose values are given in
Table 1.
At high ionic strength (/ > 0.1 m ) activity coefficients are
not simply universal functions of ionic strength, but are dependent on specific ion interactions (Felmy and Weare, 1986).
Activity coefficients are now expressed in a virial-type expansion of the form (Pitzer, 1979)= ln-y?H + 2

1.10(Ca2+)/(Mg2+)+ blwhere A and B are constants depending only on the dielectric
constant, solution density, and temperature, z is the ioniccharge, a and b are adjustable parameters fit to mean salt
activity coefficients, and / is the ionic strength. The adjustableparameters a and b for individual species are listed in Table
1. The activities of neutral species are calculated asIn-/ = a'I

1 + 10.0- 2+ 5
(Ca )°- /(Na+)

XC<•aclay

Activity Coefficients
At low to moderate ionic strength, we calculate activity
coefficients with an extended version of the Debye-Hiickel
equation (Truesdell and Jones, 1974)

4.65(Ca2+)/(Mg2+)

CECom - *Mgom ~ xKm

1 +The constants 4.65 and 1.10 represent the respective Ca-Mg
selectivity constants for organic matter (Baes and Bloom,1988) and smectite (Suarez and Zahow, 1989), while the constants 10.0 and 1.96 represent the respective Ca-Na selectivity
constants for organic matter and smectite. In the absence ofsimilar data on the Ca-K exchange selectivity, we utilize theexchange expression given above in Eq. [28] for both clay andorganic matter. Use of these concepts results in a non-constantselectivity constant for the bulk soil, with decreased Ca selectivity with increasing exchangeable Ca values (at the low rangeof exchangeable Ca values, consistent with experimental data).
We utilize this exchange model whenever the soil organic
matter content is known.

[30]

)

~v

V^E,^ora — -^Naomt Truesdell and Jones, 1974.
i Reardon and Langmuir, 1976.

In -y = —The system is defined by the following set of equations:
_

-0.5
-0.5

H,SiOS+ . . . [39]where m, is the molality, y"H is a modified Debye-Huckel
activity coefficient, which is a universal function of ionicstrength, and By and Qfc are specific coefficients for eachion-ion interaction. The specific form of this equation for
cations, anions, and neutral species is adopted from Felmyand Weare (1986).

1638

SOIL SCI. SOC. AM. J., VOL. 61, NOVEMBER-DECEMBER 1997

Table 2. Temperature-dependent thermodynamic constants.
Constant
log K COZ
log K,, H2CO3
log KBl H2CO,
log K CaCO°3
log K CaSOS
log K CaHCO,+
log K MgCO°3
log K MgSO,
log K MgHCOJ

«i

«2

a,

a,

«s

References

108.3865-356.3094-107.8871
1228.732-1.24-6 919.53
21 834.37
5 151.79-35 512.75
034 765.05-3265
0-11 132.0
93.64230.77024 470.99
2 839.319

0.01985076-0.06091964-0.03252849
0.299444-0.0036-0.31294-0.04467-0.011-0.13380.001830
0.0006058
0.01705-0.077993-40.45154
126.8339
38.92561-485.8180478.78200000071.595

669 395-1684915-563 713.9
0
0
0
0
0
0
0
0
0
0

Plummet and Busenberg (1982)
Plummer and Busenberg (1982)
Plummet and Busenberg (1982)
Plummer and Busenberg (1982)
Bell and George (1953)
Plummer and Busenberg (1982)
Reardon and Langmuir (1974)
Jacobsen (1973)
Reardon (1974, p. 85)
Robinson and Stokes (1965)
Robinson and Stokes (1965)
Stumm and Morgan (1981)
Plummer and Busenberg (1982)-1209.120

21.39
0.95
76.344-1.15083-0.76645
6.0875-171.9065

log A
log B
log K.
log A' calcite

Temperature Dependence
Most of the equilibrium constants used in this model havebeen determined as a function of temperature. Where available, the equilibrium constants are expressed as a power function of the absolute temperature:[40]

log* =

Table 2 lists the temperature-dependent thermodynamic constants used in the calculations.
The temperature dependence of the equilibrium values forwhich the constants of Eq. [40] do not exist are calculated
from the enthalpy of reaction (Table 3).

Kinetic Model of Calcite Precipitation-Dissolution
Calcite precipitation-dissolution kinetic models are almost
all based on the assumption that the reaction rate is dependent
on the surface area of the calcite. The dissolution-precipitation
model of Plummer et al. (1978) is the most comprehensivemodel for pure systems, based on dissolution studies across awide range in pH (2-7.0) and CO2 pressure (35 Pa-100 kPa).
The calcite dissolution in the absence of dissolved organiccarbon (DOC) is thus calculated with the rate equation of
Plummer et al. (1978):Rc = k,and k{ = 0.051. The dissolution-precipitation rate Rc is expressed in millimoles of calcite per square centimeter of surface area per second. The term (Us) is the H + activity at thecalcite surface. It is assumed to be (H+) at calcite saturation
where pCO2 at the surface equals pCO2 in the bulk solution.
The temperature dependence of the constants k\, k2, and k2
is expressed as+ k2 (H2CO3*) + A:3 (H2O)(Ca2+)(HCOf )k, = k( +where at 25°C, k, = 0.051, k2 = 3.45 X 10~5, k, = 1.18 X 10~7,
Table 3. Temperature dependence of equilibrium constants[43]

where the constants are given in Table 4. For conditions where
pH > 8 and pCO2 < 1.01 kPa, the Plummer et al. (1978)
equation underestimates the precipitation rate (Suarez, 1983;
Inskeep and Bloom, 1985); thus the following expression isused (Inskeep and Bloom, 1985):

R'c = -11.82[(Ca2+)(COn -[44]

with an apparent Arrhenius activation energy of 48.1 kJ mol~'for the rate constant.
The above relationships for calcite crystal growth are basedon a clean calcite surface in the absence of "surface poisons".
The inhibiting effects of DOC on calcite precipitation is wellestablished and related to surface adsorption of DOC. Inskeep
and Bloom (1986) reported data on the effect of water-solubleorganic C on calcite crystal growth. From their data, we developed the following predictive relationship for reduction in
calcite precipitation due to the presence of DOC:

[41]

(H2C03*) + k, (H20)] [42]I = 1, 2, 3

log kj =a2x2 —

r =[45]

where r is the reduction constant, x is the DOC (jxmol LT1)and «!, a2, and a3 are the nonlinear regression coefficients(0.005104,0.000426, and 0.069111, respectively). This relation,
as shown in Fig. 2, provided excellent prediction, with a correlation coefficient of 0.998. The reduction constant, r, is thenmultiplied by the Rc values calculated with either Eq. [41] or(Truesdell and Jones, 1974).
SpeciesK(289.15°K)

J mol-'
AqueousNaCOf
NaSOr
NaHCO?
KSO,H,Si04HjSiOfT
B(OH)4Solid
GypsumHydromagnesite
Nesquehonite
Sepiolitet K, association.
t Ball et al., 1987.

K.t
0.05400.12000.562
0.1413
1.175 x 10-'°2.404 x 10-22
5.75 x 10-'°Ksf
2.512 x 10-5$
1.730 x 10-37
2.393 x 10^'38
3.138 x 10--37337
-4693
-9428
37437
124502-42 400
1131-106930-28 366-

Table 4. Temperature dependence of kinetic constants of calcite(Plummer et al., 1978) and dolomite (Busenberg and Plummer, 1982).
Calcite0.198
k, (T < 298.15)
fcj (T > 298.15)

2.84
-5.86
-1.1
Dolomite

*4t

2.12
-0.07
0.53
3.16

k.

7.34

k,

t For FeCO, < 0.007.

-444
-2177
-317
-1737
-1880
-1800
-2700
-2300
-3700

1639

SUAREZ & S1MUNEK: UNSATCHEM MODEL FOR UNSATURATED TRANSPORT

1

RD = /c1(H+)05 + A:2(H2CO3*)05

y = exp(-0.005104* - 0.000426*2 - 0.0691 \xm)

+ A:3(H2O)05 - &4(HCO3~)

0.9
0.8

x Observed

0.7— Fitted

sO-6
c

|0.5
•o

"0.4
0.3

0.2
0.1

20

40

60

SO

100

Dissolved organic carbon (DOC) [>*mol/l]

120

140

160

Fig. 2. Relationship between reduction in calcite precipitation anddissolved organic carbon (DOC). Data are taken from Inskeep
and Bloom (1986).
[44] to obtain the predicted rate constant in the presence ofthe specified DOC concentration.
These and other calcite rate models all consider reaction
rates to be proportional to surface area. For simulation of
calcite dissolution in natural systems, these models may besuitable, after adjustment for the poisoning of the surface,
as discussed above. However, these rate models may not besuitable for predicting calcite precipitation rates, as the concentrations of DOC in natural environments are usually comparable to levels found by Inskeep and Bloom (1986) to completely inhibit calcite crystal growth.
Recently Lebron and Suarez (1996) developed a precipitation rate model that considers the effects of DOC on hetergeneous calcite nucleation. In the presence of clay, the rate RT isequal to the sum of the crystal growth rate and heterogeneous
nucleation rate, where

RCG = fcCGSc[(Ca2+)(COn - KS
X [0.14 - 0.11 log DOC]

[46]

and

Rm = km0.02S log[n - 2.5][3.37 X 1Q-U4DOC] [47]where kca - 64.8 mmol s~' rrT2, Sc is the calcite surface area
(m2 L-'), ArHN = 7.82 X 10~4 s"1 m~ 2 , S is the total soil surface
area (m2 L"1), DOC is expressed in mmol L~', fl = IAP/KSP,
and /?HN is the precipitation rate (mmol L~' s'1). The value
0.025 represents the active sites related to the total surface
area. This value is given for illite and may vary with mineralogy. This equation leads to calcite precipitation rates that are
independent of the calcite surface area when DOC >0.05 mM,
consistent with the experimental data of Lebron and Suarez(1996). This model considers heterogeneous calcite nucleation(on clay surfaces) with poisoning of the crystals upon formation. We consider Eq. [47] to be the most realistic for precipitation in root-zone environments while Eq. [44] and [45] aremost suitable for dissolution. The use of the Plummer et al.
(1978) equation is suitable only for very low DOC environments (<0.01 mM).

Kinetic Model of Dolomite Dissolution
The reaction rates of dolomite dissolution were calculatedwith the rate equation of Busenberg and Plummer (1982) for
sedimentary dolomites:[48]

where the temperature dependence of the constants kt, k2, k},
and fc, is expressed by Eq. [43] (Table 4). The precipitation
rate RD is expressed in millimoles of dolomite per squarecentimeter of dolomite surface area per second. According toBusenberg and Plummer (1982), the reaction rate became so
low above IAF values of 10~20 that the rate could not bemeasured within the time frame of their experiments. In theabsence of data at low levels of undersaturation, we extrapolate their results to IAP = 10~19 and assume that the reaction
rate is negligible above that value. As a result, dolomite equilibrium in the model is never achieved by dissolution, butequilibrium or supersaturated values can result from decreases
in pCO2 or by concentration of salts by evapotranspiration.
In the absence of data specific to dolomite, we utilize thecalcite reduction constants for the effect of organic matter on
dissolution rates.
Kinetic Model of Silicate Mineral Dissolution
Although silicate mineral dissolution rates are orders ofmagnitude slower than calcite dissolution, they nonethelessare extremely important for predicting base cation concentrations in noncalcareous soils and for assessing sensitivity tosoil acidification, either by natural processes or by acid rain.
Several different rate expressions have been used for feldspar
dissolution, the most successful of which are variations of the
Furrer and Stumm (1986) model:

RL = kLCl[49]

where RH and RL are the proton- and ligand-promoted rates,
kL and kH are the rate constants, C& and Cl are, respectively,
the surface concentrations of protons and ligands, and n is
the order of the reaction. This model, developed for oxides,
is not able to simulate the silicate dissolution rates in the pH
range 3 to 8 (Amrhein and Suarez, 1988). Amrhein and Suarez(1988) modified this equation by adding a rate term proportional to the uncharged surface silanol groups and substituting
the term F, which represents the sum of the proton and hydroxyl sites, for the term C^ in Eq. [49]. Combining the surface
proton and hydroxyl sites into one expression was justified by
the essentially equal effects of the two surface groups on
the dissolution rates of plagioclase feldspars. This has been
determined experimentally in the work of Chou and Wollast(1985) for albite and Amrhein and Suarez (1988) for anorthite,
representing the end members in the plagioclase feldspar series. However, other studies indicate that this equality cannot
be generalized, thus we use the modified rate equation
4= a(r)[50]

where a, a', £>,, and b2 are the rate coefficients for the proton,hydroxyl, neutral, and surface-ligand sites and F, F', SOH,
and S-L are the surface concentrations of proton, hydroxyl,
neutral, and surface-ligand sites. Detailed rate data for use
with Eq. [50] are not available for all silicates; however, inthe pH range of 5 to 9, this equation can often be simplified
to#, = ^(SOH) + £2(S-L)
8

[51]
6

where fc, and b2 are equal to 2.09 X 10~ and 4.73 X 10~ s'1,
respectively, for anorthite, and SOH, the total number ofsurface sites, is taken as 0.12 mmol m~ 2 (Amrhein and Suarez,

1640

SOIL SCI. SOC. AM. J., VOL. 61, NOVEMBER-DECEMBER 1997

1988) when R, is expressed in mmoles per square meter per
second.
Using the rate model expressed in Eq. [51], and in that
restricted pH range, we consider the following reactions foranorthite, labradorite, albite, K feldspar, biotite, and hornblende:

CaAl2Si2O8 + H2CO3 + H2O -*Ca2++ 2HCO3~ + kaolinite[52]

4Cao.5Nao.5AlL5Si2.5Og + 6H2CO3 + 11H2O -»Values for these silicate minerals are taken from the studies
of Swoboda-Colberg and Drever (1993) and Suarez and Wood
(1996), based on long-term weathering studies in experimentsthat minimized the effects of grinding and pretreatment artifacts.
An alternative approach is the rate model of Sverdrup andWarfvinge (1988):
frTJ

r\

Morel

"

[58]
where k is the reaction rate for the different processes, andn, x, y, u, and m are the reaction orders, all determined experimentally. We utilize the field values given by Sverdrup andWarfvinge (1988).
Adsorption
Transport of minor element anions such as B, As, and Se
is often affected by adsorption processes. Various modelingapproaches have been used in the past, ranging from retardation factors to the use of Langmuir isotherms. These models
generally require the input of adsorption isotherm values specific to a given pH. In this model, we utilize the constantcapacitance model, since it is able to represent the effect ofpH on adsorption affinity, and our model has the capacity to
predict soil solution pH. The model contains the following
assumptions: adsorption is a ligand-exchange mechanism, allsurface complexes are inner sphere, and no surface complexes
are formed with other salts in solution. The relation between
the surface charge and the surface potential is given by
o- =

[60]

K=

-

=

2Ca2+ + 2Na+ + 6HCO3- + Skaolinite + 4H4SiO4[53]
2NaAlSi3O8 + 2H2CO3 + 9H2O ->2Na + + 2HCO3- + kaolinite + 4H4SiO4[54]
2KAlSi3O8 + 2H2CO3 + 9H2O ->2K+ + 2HCO3- + kaolinite + 4H4SiO4 [55]2KMg3AlSi3O10(OH)2 + 14H2CO3 + H2O -»2K++ 6Mg2+ + 14HCO3- + kaolinite + 4H4SiO4 [56]Ca2Mg4Al2Si7O22(OH)2 + 12H2CO3 + 5H2O —2Ca2+ + 4 Mg2+ + 12HC03- + kaolinite + 5H4SiO4[57]

FTLJ + ln

is expressed in molc L '. The intrinsic conditional equilibrium
constants corresponding to B adsorption are (Goldberg and
Glaubig, 1985)

B[SH2B03][SOH][H3B03][62]

The mass balance equations for the surface functional groups
and B are given by[SOH]T = [SOH] + [SOH2+] + [SO'] + [SH2BO3][63]
BT = H3BO3 + SH2BO3[64]

and the charge balance equation for the surface is defined by
CT = [SOH2+] - [SO'][65]

In the absence of soil-specific data, we utilize the average soil
constants of 9.3, -10.6, and 5.5 for the values of log K+, log
K-, and log KB, respectively (Goldberg, 1993). We utilize thevalue of 1.06 F m"2 for the capacitance (C). Using the dataof Goldberg and Glaubig (1985), we developed the following
expression relating the experimentally determined adsorption
site density to S, the soil surface area:[SOH]T = 2.53 X ICT7 + 4.61 X lO'9 S[66]

where S is expressed in square meters per gram and adsoption
density is in moles per gram. As shown in Fig. 3, this expression
provides a satisfactory prediction of B adsorption density as
a function of total surface area. In contrast to direct determination of adsoprtion density, which requires time-consuming
adsorption experiments, surface area can be determined relatively rapidly or estimated from mineralogy and clay content.
The relation given by Eq. [66] had an r value of 0.96 (n =14). By inputting a soil surface area and calculating suspension
density a from water content and bulk density, the system is
now defined by seven equations and seven unknowns.
The numerical solution of the chemical system is madeusing an iterative approach similiar to that utilized in WATEQ(Truesdell and Jones, 1974), except with the addition of solidand sorbed-phase reactions. Initial calculation is made of ionicstrength, temperature-dependent constants, and activity coefficients. Calculation is then made of pH from charge balanceand partioning of alkalinity into HCO3 and CO3 by solving
Eq. [Al] to [A4] and equilibrating the solution composition
with the solid phases. The iteration criterion in this section is
the second dissociation of H2CO3. After passing this criterion,
the major ion concentrations are redistributed into the ion
pairs and complexes. Each component is checked for massbalance (i.e., Eq. [12]). After iteration, the sorbed phases
are now equilibrated with the solution, and mass and chargebalance checked. The process is repeated until specified convergence is achieved (typically set at 0.1%). Overall convergence criteria include charge balance and change in ionicstrength.

[59]

where C is the capacitance density (F m~ 2 ), S is the specificsurface area (m2 g"1), a is the suspension density (g L~'), F
is the Faraday constant (C molc~'), ^ is potential (V), and a

RESULTS AND DISCUSSION
Unique features of the model are demonstrated in

the following simulations. In the first simulation, weSUAREZ & SIMUNEK: UNSATCHEM MODEL FOR UNSATURATED TRANSPORT

1641

y=2.53.10-7+4.61
1.4E-620

1.2E-6W

1 .OE-6

I

O

•o

40

8.0E-7cO•

6.0E-7W

4.0E-760.Q80
200-365 d

2.0E-7

20

40

60

80

100

120

140

160

180

200

Soil surface area (m2 g~ 1 )

100
O.I

Fig. 3. Relationship between B absorption density and soil surface
area for 14 western U.S. arid-zone soil samples. Data taken fromGoldberg and Glaubig (1985).

examine differences in solution composition and location of precipitation when simulating calcite precipitation using a kinetic as compared with an equilibrium
model. We specify irrigation with a Colorado River derived well water, which is already calcite supersaturated(drainage well no. 15 from the Wellton-Mohawk Irrigation District in Arizona, as analyzed by Suarez [1977a]).
The soil CO2 distribution is taken as atmospheric at thesoil surface, increasing linearly to 2 kPa (2%) at thebottom of the root zone. In this example, the water
uptake distribution by plants is also linearly scaled fromthe maximum at the soil surface to zero at the bottomof the soil root zone (100 cm). The initial volumetricwater content is 0.143 m3 m~3. We assume a constantirrigation of 1 cm d~' and a plant water uptake of 0.9
cm d"1. In this example, we neglect ion exchange inorder to evaluate differences due to the kinetic andequilibrium assumptions. The unsaturated hydraulicconductivity-water content relationship was calculatedusing the equation of van Genuchten (1980) from the

O.Z

0.3

0.4

Water Content (9)Fig. 4. Soil water content profiles with depth for 365-d simulation ofirrigation and plant water uptake using Wellton Mohawk drainage
water as the irrigation source.

water content-pressure head data given in Hillel andvan Bavel (1976).
Figure 4 shows that the water content in the soil
increased with time and was close to steady state after
100 d. The decrease in water content with depth is due
to plant root extraction of water. The Ca concentration
profile using the equilibrium model, shown in Fig. 5a,increases steadily with time, reaching steady state after
100 d. After 100 d, the concentration ranges from 6
mmolc L~' at the surface to 32 mmolc L"1 at the 60-cm
depth. Using the kinetic model with an assumed calcitesurface area of 1 X 10~7 m2 m'3 of soil, Ca concentrations
are almost doubled near the surface and similar to theequilibrium case at the bottom of the profile. The Cavalues are roughly similar at the bottom of the soil
profile because of gypsum precipitation at the lower
depths, discussed below. The simulated output for the
kinetic model is dependent on the surface area value

o

H

Q.
o>

O)

Q

Q

80 -

100

Ca Concentration [mmolc

Ca Concentration [mmolc L"1]

Fig. 5. Calcium concentration profiles with depth for 365-d simulation of irrigation and plant water uptake using Wellton Mohawk drainage
water and (a) calcite equilibrium condition and (b) kinetic mode) with calcite surface area taken as 1 x 10 7 m2 m ' of soil.

1642

SOIL SCI. SOC. AM. J., VOL. 61, NOVEMBER-DECEMBER 1997

I
I

Q

80 -

IOO

Alkalinity [mmolc L" ]Fig. 6. Alkalinity concentration profiles with depth for 365-d simulation of irrigation and plant water uptake using Wellton Mohawk
drainage water and calcite equilibrium condition compared withkinetic model with calcite surface area taken as 1 x 10 7 nr m 'of soil.

and DOC content of the soil water. In this case, we did
not include DOC.
The alkalinity values shown in Fig. 6 indicate large
differences between the two simulations, with the kinetic model again predicting much larger concentrations, especially in the shallow depths but still twiceas large at the bottom of the profile. These alkalinity
differences are related to the calcite supersaturation
throughout the soil profile, as shown in Fig. 7, whichshows the pIAP values for calcite as a function of depth
and time. The final steady-state pIAP values range
around 7.0 in the near-surface environment to 8.1 at
the 100-cm depth, similar to those reported by Suarez(1985). In contrast, the equilibrium model results in

20—>
Equilibrium

,o
.c
S.

CD

Q

40

60

80

100

6.5

7.5

8.5

plAP,(Ca2+)(CO,Fig. 7. Calcium carbonate pIAP (ion activity product) values withdepth for 365-d simulation of irrigation and plant water uptake
using Wellton Mohawk drainage water. Line at pIAP = 8.48 represents value for equilibrium condition and curved lines represent
values obtained using the kinetic model with a calcite surface areaof 1 x 10~7 m2 m 3 of soil.

calcite equilibrium (pIAP = 8.47) at all times and locations in the profile.
Due to the initial equilibrium assumption, the kinetic
simulation took longer to reach steady state, as evident
by the pIAP-time relations in Fig. 7, as well as thecomparison of alkalinity-time relations in Fig. 6. In addition to differences in solution composition, the models
also give differences in the carbonate distribution withdepth, as shown in Fig. 8a and 8b. The equilibrium
model predicts a large quantity of calcite precipitationin the first depth followed by a very large quantity ofcalcite dissolution in the shallow depths, followed byslight precipitation at depth. This unrealistic result demonstrates the limitations of an equilibrium model for
simulation of calcite precipitation. The result can bepartially avoided by the equally unrealistic forcing of
the input solution to calcite saturation before initiating
the program. In contrast to the equilibrium model, thekinetic model simulates a condition of decreasing precipitation with depth. The predicted distribution depends on various factors including the calcite saturation
status of the input water, amount of water applied, CO2
distribution with depth, plant water uptake and distribution with depth, as well as potential gypsum precipitation, which in this case limits calcite precipitation at
depth. Quantity and depth of calcite precipitation in soil
profiles has been used in various simplified equilibrium
models for reconstruction of paleoclimate as well as forpedogenesis. While simplified models may be needed
for regional-scale climatic simulations, nonetheless, adetailed mechanistic model provides the opportunity to
examine the sensitivity of the results to input assumptions. In addition this example illustrates the importanceof kinetic factors on calcite distribution — equilibrium
models will not yield realistic distributions.
Gypsum precipitation using the calcite kinetic model
is simulated as occurring in the 60- to 100-cm zone, witha maximum at 67 cm. In contrast, with the equilibrium
model, the location of gypsum accumulation is shifted
downward, due to the smaller Ca and alkalinity concentrations for the equilibrium than for the kinetic model.
The second example considers both the interactionof transient water flow and^ water content, as well as
the dynamic CO2 submodel (Simunek and Suarez, 1993),
with the equilibrium and kinetic carbonate subroutines.
In these simulations we considered intermittent irrigation of 10 cm, every 10 d, of Colorado River water(where Ca = 2.63, Mg = 1.03, Na = 2.55, Cl = 1.94,
alkalinity = 2.33, and SO4 = 2.0 mmolc LT')- Although
the model can calculate CO2 production as a functionof environmental stress, for the purposes of comparing
the different simulations we fixed the CO2 production
at 0.007 m3 m~ 2 , which is a realistic field value (Suarez
and Simunek, 1993). Other conditions are those given
for the previous example.
After 100 d, a steady-state cycle was reached in thatthe water content and solution composition cycled between irrigations. Shown in Fig. 9 is the water content
distribution with depth for a 10-d interval between irrigations after 300 d of simulation. After irrigation, thewater content increases in the shallow soil zone with

SUAREZ & SIMUNEK: UNSATCHEM MODEL FOR UNSATURATED TRANSPORT

100
-40

-20

O

20

1643

40

1

Calcite [mmol kg " soil]
Calcite [mmol kg -1 soil]
Fig. 8. Quantity of calcite precipitated (+) or dissolved (-) with depth for 365-d simulation of irrigation and plant water uptake using Wellton
Mohawk drainage water and (a) calcite equilibrium condition and (b) kinetic model with calcite surface area taken as 1 x 10 7 nr m~ 3 of soil.

subsequent redistribution and water uptake by plants.
Ten days after an irrigation, the water content is backto the initial conditions immediately before the next
irrigation. Using the dynamic CO2 production transportmodel, we simulate changing CO2 concentrations withdepth and time. As shown in Fig. 10, the CO2 concentrations are greatest 2 d after an irrigation when the watercontent near the surface is still high. High water content
resticts gas diffusion, and after sufficient time restricteddiffusion results in elevated CO2 concentrations. Upondrainage, the soil CO2 distribution quickly reverts backto the pattern shown on Days 4 through 10, with increasing CO2 concentration with depth to the bottom of theroot zone. These dynamic distributions with space andtime in turn produce effects on the soil solution composition.
The simulations in Fig. 11 show that the Ca concentration is relatively constant during the 10-d cycle using

the equilibrium model at fixed CO2. Increasing concentration with depth due to root water uptake is moderated by calcite precipitation. In contrast, the Ca concentrations changed with time by more than a factor of two
in the upper portion of the profile when the kinetic
model was utilized. The kinetic model also predicts
higher concentrations throughout the profile. At shallow depths, this is due to both calcite supersaturationand fluctuations in CO2; in the lower portions of the
profile it is primarily due to calcite supersaturation during precipitation. The corresponding alkalinity concentrations for the same two simulations are shown in Fig.
12. Similar to the Ca profiles, the alkalinity concentrations and changes with time are greater using the kinetic
model than the equilibrium model. The alkalinity concentrations with time in the equilibrium simulation

O

E
.o,

Q.
Q)

.c

Q

Q.
0)

Q

80 -

100

0.10

0.20

0.30

0.40

Water Content (8)Fig. 9. Soil water content profiles with depth for a 10-d irrigation
cycle after 1-yr simulation of irrigation and plant water uptake
using Colorado River water.

co, [kPa]r

Fig. 10. Soil CO, partial pressure profiles with depth for a 10-d irrigation cycle after 1-yr simulation of irrigation and plant water uptake
using Colorado River water. Simulation uses the CO. productiontransport submodel.

1644

SOIL SCI. SOC. AM. J., VOL. 61, NOVEMBER-DECEMBER 1997

0,10 d0,10Kinetic

20

20

0,0
d

40

§

Equilibrium
<D

60

80

60

80

100
8

100

I2

10

Ca Concentration, [mmol c L' ]Fig. 11. Calcium concentration profiles with depth for a 10-d irrigation
cycle after 1-yr simulation of irrigation and plant water uptake
using Colorado River water (a) calcite equilibrium condition and(b) kinetic model coupled to CO, production-transport submodel.

model change less than the corresponding Ca changes
shown in Fig. 11, because Ca concentrations weregreater than alkalinity concentrations in the irrigationwater. This ion imbalance causes a buffering effect on
the minor ion and allows larger fluctuation of the dominant ion.
The differences in Ca and alkalinity concentration
between the equilibrium and kinetic models are consistent with the pIAP changes presented in Fig. 13 for thesame 10-d cycle. Changes in pH with depth and time
are also quite variable, as shown in Fig. 14 for the kinetic
model. This range in pH has important implications forplant nutrition as well as chemical transport.

CONCLUSION
The model presented attempts a realistic representation of the dynamic processes operating on the soil solu-.8 \0,IO
x

*•

20

d

„.

9.5

8.5

7.5

10

pIAP, (Ca^ (CO 32-\
Fig. 13. Calcium carbonate pIAP (ion activity product) values withdepth for a 10-d irrigation cycle after 1-yr simulation of irrigationand plant water uptake using Colorado River water and kinetic
model coupled to CO, production-transport submodel. The value
pIAP = 8.48 represents calcite equilibrium condition.

tion composition of arid land and nonacidic soils. It
contains several unique features such as the predictionof CO2 pressure, pH, inclusion of the constant capacitance model for pH-dependent B adsorption, consideration of the effect of solution composition on hydraulicproperties, and inclusion of kinetic rather than equilibrium expressions for reactions where equilibrium cannot
be assumed. Consideration of both calcite kinetics andCO2 production and transport was demonstrated to haveimportant effects on solution composition as well as on
the quantity and location of inorganic C precipitationdissolution in the soil profile. The dynamic nature of
the soil solution composition also results in differencesin the changes experienced by different ions, and differences in changes depending on depth. Most soil solution
sampling strategies are not designed to capture thisvariation.
Some aspects of the model can be considered concep-.

Kinetic

20
Equilibrium

o 40

§

40

60

60

80

80

100

100

6.5

Alkalinity, [mmol c L
Fig. 12. Alkalinity concentration profiles with depth for a 10-d irrigation cycle after 1-yr simulation of irrigation and plant water uptake
using Colorado River water (a) calcite equilibrium condition and(b) kinetic model coupled to CO2 production-transport submodel.

7.0

7.5

8.0

8.5

PH

Fig. 14. Soil water pH values with depth for a 10-d irrigation cycle
after 1-yr simulation of irrigation and plant water uptake using
Colorado River water and kinetic model coupled to CO, production-transport submodel.

SUAREZ & SIMUNEK: UNSATCHEM MODEL FOR UNSATURATED TRANSPORT

tual in that individual parts have been determined experimentally but they have not yet been evaluatedwithin an integrated system. In addition to needed evaluation of the model with field studies, it is our objective
to incorporate into the model additional solid-phasereactions including evaporite minerals.
APPENDIX
The chemical system includes nine major aqueous components [Ca2+, Mg2+, Na + , K + , CO2.-, SO5r, CT, B(OH)3, andH4SiO4], 13 complexed species [CaCO§, CaHCOf, CaSOS,
MgCO§, MgHCOr, MgSOS, NaCOr, NaHCO?, NaSO4~,
KSO4-, B(OH)4-, H3SiO4-, and H2SiO^], six precipitated species (calcite, gypsum, dolomite, nesquehonite, hydromagnesite, and sepiolite), five sorbed species (Ca, Mg, Na, K, andSH2BO3), and five additional species (pCO2, H2CO3*,
HCOf, H + , and OH~). Altogether we consider 38 species,
and therefore we define 38 independent equations to solve
this system.
Carbon Dioxide-Water System
The dissociation reaction of water is written as

H,O ^ H+ + OH'_ (H + )(OH-)

KM =

(H20)

[Al]

where A"w is the dissociation constant for water and brackets( ) denote ion activity.
The solubility of CO2(g) in water is given by
H2CO3* ^ CO2(g) + H2O
=
C 2

°

(H2CO?) = (H 2 COf)
(C02(g))(H20)
pCo2(H20)

f
L

,
J

where the activity of CO2(g) is expressed in terms of the partialpressure pCO2, KCo2 is Henry's law constant, and H2CO3* represents both CO2 aqueous and H2CO3. Protolysis reactions of
dissolved CO2 are written as
H2CO? ^ H + + HC03HC03-

cor-H+ +. _ _ ,Ka] = (H(^Q("+) [A3]x».^Qep-„,(HCO-)[A4]

where Ka] and K,2 are the first and the second dissociation
constant of H2CO3.

1645

1646

SOIL SCI. SOC. AM. J., VOL. 61, NOVEMBER-DECEMBER 1997