The UNSATCHEM Modulefor HYDRUS (2D/3D)
Simulating Two- and Three-Dimensional Movement ofand Reactions Between Major Ions in Soils
Version 2

Jirka Šimůnek1, Miroslav Šejna2, and M. Th. van Genuchten3
May 2021
1

Department of Environmental Sciences
University of California Riverside
Riverside, CA, 92521, USA
2

PC Progress, Ltd.
Prague, Czech Republic
3

Department of Mechanical Engineering
Federal University of Rio de Janeiro
Rio de Janeiro, Brazil© 2021 J. Šimůnek and M. Šejna. All rights reserved.

2

Table of Contents
Table of Contents................................................................................................................................. 3
List of Figures ...................................................................................................................................... 5
List of Tables ........................................................................................................................................ 7
Abstract ................................................................................................................................................ 9
1. Introduction ............................................................................................................................11
2. Multicomponent Solute Transport .......................................................................................... 13
3. Carbonate Chemistry ............................................................................................................15
3.1.

Principles .....................................................................................................................15

3.2.

Mass and Charge Balance Equations ..........................................................................16

3.3.

CO2 - H2O System ........................................................................................................17

3.4.

Complexation Reactions ..............................................................................................17

3.5.

Cation Exchange and Selectivity .................................................................................18

3.6.

Precipitation-Dissolution Reactions ............................................................................19

3.7.

Kinetic Model for Calcite Precipitation-Dissolution...................................................20

3.8.

Kinetic Model of Dolomite Dissolution .......................................................................21

3.9.

Silica Concentration in Soil Solution...........................................................................22

3.10.

Activity Coefficients .....................................................................................................22
3.10.1. Extended Debye-Hückel Expression ................................................................22
3.10.2. Pitzer Expressions ............................................................................................23

3.11. Temperature Effects......................................................................................................24
3.12. Osmotic Coefficient and Osmotic Pressure ..................................................................24
3.13. Effect of Solution Composition on Hydraulic Conductivity .........................................24
4. Preprocessing..........................................................................................................................29
4.1.

Main Processes ............................................................................................................29

4.2.

Solute Transport - General Information ......................................................................29

4.3.

Solution Compositions for the UNSATCHEM Module ..................................................30

4.4.

Solute Reaction Parameters .........................................................................................31

4.5.

Chemical Parameters for the UNSATCHEM Module ...................................................32

4.6.

Initial Conditions ..........................................................................................................33

5. Post-processing .......................................................................................................................35
5.1.

Results – Graphical Display ........................................................................................35

3

5.2.

Results – Other Information.........................................................................................36

6. Input and Output Files ..........................................................................................................37
6.1.

Pitzer's Coefficients .....................................................................................................37

6.2.

The Selector.in Input File ............................................................................................37

6.3.

UNSATCHEM Output Files .........................................................................................41

7. Example Problems .................................................................................................................43
7.1.

Example 1 - Column Infiltration....................................................................................43

7.2.

Example 2 - Two-Dimensional Furrow Irrigation Problem ..........................................47

References .....................................................................................................................................53

4

List of Figures
Figure 4.1.

The Main Processes dialog window for the UNSATCHEM module. ....................29

Figure 4.2.

The Solute Transport dialog window for the UNSATCHEM module. ...................30

Figure 4.3.

The Solute Reaction Parameters dialog window for the UNSATCHEM module. .32

Figure 4.4.

The Solution Compositions dialog window for the UNSATCHEM module. .........31

Figure 4.5.

The Chemical Parameters dialog window for the UNSATCHEM module. ...........33

Figure 4.5.

The "Initial Conditions" part of the Data Tab of the Navigator Bar for the
UNSATCHEM module. ..........................................................................................33

Figure 5.1.

The "Results - Graphical Display" part of the Data Tab of the Navigator Bar for
the UNSATCHEM module. ....................................................................................35

Figure 5.2.

The Chemical Mass Balance Information dialog window for the UNSATCHEM
module. ....................................................................................................................36

Figure 7. 1.

Water content profiles at various times for example 1. ............................................44

Figure 7. 2.

Tracer concentration profiles at various times for example 1. ..................................44

Figure 7. 3.

Calcium concentration profiles at various times for equilibrium (left) and b) kinetic
calcite precipitation-dissolution for example 1. ........................................................45

Figure 7. 4.

Alkalinity profiles at various times for equilibrium (left) and kinetic (right) calcite
precipitation-dissolution for example 1. ...................................................................45

Figure 7. 5.

Calcite profiles at various times for equilibrium (left) and kinetic (right) calcite
precipitation-dissolution for example 1. ...................................................................46

Figure 7. 6.

Schematic representation and finite element mesh of the flow domain for the furrow
irrigation system for example 2. ...............................................................................48

Figure 7. 7.

Pressure head (cm) profiles at times: a) 0.1, b) 0.5, c) 1, and d) 2 days for example 2.
.................................................................................................................................49

Figure 7. 8.

Tracer concentration (-) profiles at times: a) 0.1, b) 1, c) 3, and d) 5 days for example
2. ..............................................................................................................................50

Figure 7. 9.

Sodium concentration (mmolcL-1) profiles at times: a) 0.5, b) 1, c) 3, and d) 5 daysfor example 2. ..........................................................................................................51

Figure 7. 10. Exchangeable concentrations of sodium (mmolckg-1) profiles at times: a) 0.5, b) 1,
c) 3, and d) 5 days for example 2. ...........................................................................525

6

List of Tables
Table 3.1. Commands in the Project Manager. .............................................................................15
Table 5.1. Additional UNSATCHEM variables displayed in the View Window of the Resultstab (Results - Graphical Display).. ...............................................................................35
Table 6.1. Major ion chemistry information in the Selector.in input file. ....................................38
Table 6.2. EQUIL.OUT - chemical information. ..........................................................................427

8

AbstractŠimůnek, J., M. Šejna, and M. Th. van Genuchten, The UNSATCHEM Module for HYDRUS(2D/3D) Simulating Two- and Three-Dimensional Movement of and Reactions Between Major
Ions in Soils, Version 2.0, PC Progress, Prague, Czech Republic, 54 pp., 2021.
This report documents the geochemical UNSATCHEM module [Šimůnek and Suarez, 1994;
Šimůnek et al., 1996] that has been implemented into the two- and three-dimensional
computational module of the HYDRUS (2D/3D) software package. The geochemical
UNSATCHEM module simulates the transport of major ions in variably-saturated porous media,
including major ion equilibrium and kinetic non-equilibrium chemistry. The resulting code isintended for predictions of major ion chemistry, along with water and solute fluxes in soils
during transient flow. The major variables of the chemical system in UNSATCHEM are Ca, Mg,
Na, K, SO4, Cl, NO3, H4SiO4, alkalinity, and CO2. The model accounts for various equilibrium
chemical reactions between these components, such as complexation, cation exchange, andprecipitation-dissolution. For the precipitation-dissolution of calcite and dissolution of dolomite,either equilibrium or multicomponent kinetic expressions can be used, which includes bothforward and back reactions. Other dissolution-precipitation reactions considered include gypsum,hydromagnesite, nesquehonite, and sepiolite. Since the ionic strength of soil solutions can varyconsiderably in time and space and often reach high values, both the modified Debye-Hückel and
the Pitzer expressions are incorporated into the model to calculate single ion activities. The effectof solution chemistry on the hydraulic conductivity is also considered. Water flow and heattransport modules are similar (almost identical) as in regular HYDRUS. Applications of the
UNSATCHEM module are demonstrated in several examples.
This report serves as both a Technical and User Manual and is a reference document of theGraphical User Interface of the UNSATCHEM related parts of the HYDRUS software package.

9

DISCLAIMER
This report documents the UNSATCHEM module of HYDRUS (2D/3D). The UNSATCHEM
module was developed as a supplemental module of the HYDRUS software package to model thetransport and reactions of major ions in the soil. The software has been verified against selected
test cases. However, no warranty is given that the program is completely error-free. If you do
encounter problems with the code, find errors, or have suggestions for improvement, please contact
one of the authors at
Jirka Šimůnek
Tel/Fax: 1-951-827-7854
Email: jiri.simunek@ucr.edu

10

1. Introduction
The main purpose of this report is to document the UNSATCHEM module [Šimůnek and Suarez,
1994; Šimůnek et al., 1996] of the HYDRUS (2D/3D) software package [Šimůnek et al., 2011;Šejna et al., 2011] for simulating two- and three-dimensional multicomponent transport andreactions of major ions. The two-dimensional variably-saturated water flow and heat movement aredescribed in the HYDRUS (2D/3D) documentation and will not be repeated here. The
UNSATCHEM module is fully supported by the HYDRUS (2D/3D) graphical user interface [Šejna
et al., 2011].
The resulting coupled program numerically solves the Richards equation for saturated-unsaturated
water flow and convection-dispersion type equations for heat and solute transport. The water flow
equation incorporates a sink term to account for water uptake by plant roots. The heat transport
equation considers movement by conduction as well as convection with flowing water. The major
variables of the chemical system are Ca, Mg, Na, K, SO4, Cl, NO3, H4SiO4, alkalinity, and CO2. The
model accounts for equilibrium chemical reactions between these components, such ascomplexation, cation exchange, and precipitation-dissolution. For the precipitation-dissolution of
calcite and dissolution of dolomite, either equilibrium or multicomponent kinetic expressions areused, which include both forward and back reactions. Other dissolution-precipitation reactions
considered include gypsum, hydromagnesite, nesquehonite, and sepiolite. Since the ionic strength ofsoil solutions can vary considerably with time and space and often reach high values, both modified
Debye-Hückel and Pitzer expressions were incorporated into the model to calculate single ion
activities.
The UNSATCHEM module may be used to analyze water and solute movement in unsaturated,partially saturated, or fully saturated porous media. UNSATCHEM can handle flow domains
delineated by irregular boundaries. The flow region itself may be composed of nonuniform soils
having an arbitrary degree of local anisotropy. Flow and transport can occur in the vertical plane,
the horizontal plane, or in a three-dimensional region exhibiting radial symmetry about a verticalaxis. The water flow part of the model considers prescribed head and flux boundaries, as well asboundaries controlled by atmospheric conditions.
The governing flow and transport equations are solved numerically using standard Galerkin-typelinear finite element schemes. Applications of the UNSATCHEM module are demonstrated on
several examples.

11

12

2. Multicomponent Solute Transport
The partial differential equation governing two- and three-dimensional advective-dispersive
chemical transport under transient water flow conditions in a partially saturated porous medium istaken as∂θ ck
∂ck  ∂qi ck∂ck∂ c k∂ 
k 1, 2,..., N c
=
+ρ
+ρ
ij
 θ D=
∂t
∂t
∂t
∂xi 
∂x  ∂xi(2.1)

where ck is the total dissolved concentration of the aqueous component k [ML-3], c k is the totalsorbed concentration of the aqueous component k [MM-1], c k is the total concentration of aqueouscomponent k in the minerals, which can precipitate or dissolve [MM-1], ρ is the bulk density of themedium [ML-3], Dij is the dispersion coefficient tensor [L2T-1], qi is the volumetric flux [LT-1], andNc is the number of aqueous components. Solute uptake by plant roots is not considered in (2.1).
The second and third terms on the left side of eq. (2.1) are zero for components that do not undergo
ion exchange or precipitation/dissolution.
Dispersion coefficients and initial and boundary conditions are similar as for solute transport in thestandard module of HYDRUS (2D/3D).

13

14

3. Carbonate Chemistry
3.1. Principles
The carbonate chemistry module was adopted from the UNSATCHEM software package[Šimůnek et al., 1996], and thus their description closely mirrors material in the UNSATCHEM
manual. Additional details can be found in the original manual [Šimůnek et al., 1996].
When using the ion-association model (and Debye-Hückel activity coefficient calculations), we
assume that the chemical system for predicting major ion solute chemistry of the unsaturated
zone includes 37 chemical species. These species are divided into six groups, as listed in Table
3.1. They include 7 primary dissolved species (calcium, magnesium, sodium, potassium, sulfate,chloride, and nitrate), 10 complex aqueous species, 6 possible solid phases (calcite, gypsum,nesquehonite, hydromagnesite, sepiolite and dolomite), 4 surface species, 7 species which form
the CO2-H2O system, and 3 silica species. The species from the last two groups could have also
been included in other groups (e.g., CO32-, H2SiO42-, and H+ could be included in the first group).
Their consideration into separate groups is mainly due to their different treatment compared tothe other species. For example, complex species of these groups are also considered at high ionicstrength when the Pitzer equations are used to calculate activity coefficients, while the species ofthe second group are in that case dropped from the system, as discussed later. One of the solid
phases (dolomite) is not included in the equilibrium system since its dissolution is always treatedkinetically. Also, the exclusion of calcite from the equilibrium system is optional since its
precipitation-dissolution can be treated as a kinetic process. As a result, either 36 or 35
independent equations are needed to solve this system. In the following sections, we present thisset of equations, while the solution method is discussed in Šimůnek et al. [1996].
Table 3.1. Chemical species considered in the carbonate chemistry module.
1
2

Aqueous components
Complexed species

7
10

3

Precipitated species

6

4
5
6

Sorbed species
CO2-H2O species
Silica species

473Ca2+, Mg2+, Na+, K+, SO42-, Cl-, NO3CaCO3o, CaHCO3+, CaSO4o, MgCO3o, MgHCO3+,
MgSO4o, NaCO3-, NaHCO3o, NaSO4-, KSO4CaCO3, CaSO4⋅ 2H2O, MgCO3⋅ 3H2O,
Mg5(CO3)4(OH)2⋅ 4H2O, Mg2Si3O7.5(OH) ⋅ 3H2O,
CaMg(CO3)2
Ca,M g,N a,K
PCO2, H2CO3*, CO32-, HCO3-, H+, OH-, H2O
H4SiO4, H3SiO4-, H2SiO42-

15

3.2. Mass and Charge Balance Equations
Seven mass balance equations for the primary species in the 1st group and one for the silica
species in the 6th group of Table 3.1 are defined:
Ca T =[Ca 2+ ] + [CaSO04 ] + [CaCO30 ] + [CaHCO3+ ]Mg T =[Mg 2+ ] + [MgSO04 ] + [MgCO30 ] + [MgHCO3+ ]Na T =[Na + ] + [NaSO-4 ] + [NaCO3- ] + [NaHCO30 ]K T [K + ] + [KSO-4 ]=SO 4T =[SO 24− ] + [CaSO04 ] + MgSO04 ] + [NaSO -4 ] + [KSO -4 ](3.1)

ClT = [Cl- ]NO3T = [NO3- ]SiO 4T = [H 4SiO 4 ] + [H 3SiO-4 ] + [H 2SiO 2-4 ]in which variables with subscript T represent the total analytical concentration in solution of thatparticular species, and where brackets refer to molalities (mol kg-1). Two mass balance equations
for the total analytical concentration of carbonate and bicarbonate are defined as follows

CO3T =[CO32− ] + [CaCO30 ] + MgCO30 ] + [NaCO3- ]HCO3T =
[HCO3− ] + [CaHCO3+ ] + MgHCO3+ ] + [NaHCO30 ](3.2)

which are used to calculate inorganic alkalinity, Alk (molckg-1):
Alk = 2CO3T +HCO3T + [OH - ] − [H + ](3.3)

Most chemical and multicomponent transport models consider the total inorganic carbon to be aconservative property [e.g., Westal et al., 1986; Liu and Narasimhan, 1989; Yeh and Tripathi,
1991]. However, this approach can be used only for closed systems. In a soil profile withfluctuating CO2 concentrations, the approach is inappropriate, and the use of alkalinity as aconservative property is preferable.
In addition to the mass balance equations, the overall charge balance equation for the solution is
given as2[Ca 2+ ] + 2[Mg 2+ ] + [Na + ] + [K + ] + [CaHCO3+ ] + [MgHCO3+ ] + [H + ] − 2[CO32− ]−−[HCO3− ] − 2[SO 24− ] − [Cl- ] − [NO3- ] − [OH - ] − [NaCO3- ] − [NaSO-4 ] − [KSO-4 ]= 0

16

(3.4)

3.3. CO2 - H2O System
The activities of the species present in solution at equilibrium are related by the mass-action
equations. The dissociation of water is written as follows
H 2 O ⇔ H + OH+

(H + ) (OH - )Kw =(H 2 O)-(3.5)

where KW is the dissociation constant for water [-], while the parentheses denote ion activities.
The methods for calculating ion activities will be discussed later.
The solubility of CO2(g) in water is described with Henry's Law:

H 2 CO*3 ⇔ CO 2(g) + H 2 O

K CO2 =

(H 2 CO*3 )PCO2 (H 2 O)

(3.6)

where the activity of CO2(g) is expressed in terms of the partial pressure PCO2 (atm), KCO2 isHenry's law constant, and H2CO3* represents both aqueous CO2 and H2CO3.
Protolysis reactions of dissolved CO2 are written as

*
3

+

−
3

+

H 2 CO ⇔ H + HCO
HCO ⇔ H + CO

−
3

2−
3

(H + )(HCO3− )K a1 =(H 2 CO*3 )(3.7)(H + )(CO32− )Ka2 =(HCO3− )(3.8)

where Ka1 and Ka2 are the first and the second dissociation constants of carbonic acid [-],
respectively.

3.4. Complexation Reactions
Each complexation reaction for the species in the second group of Table 3.1 and for the silica
species can be represented by a mass action law:K6 =(Mg 2+ ) (HCO3- )(MgHCO3+ )(3.10)K5 =(Mg 2+ ) (CO32- )(MgCO30 )(Mg 2+ ) (SO 2-4 )(MgSO04 )

(3.9)

K4 =

(Ca 2+ ) (HCO3- )(CaHCO3+ )

K2 =

K3 =

(Ca 2+ ) (SO 24− )(CaSO04 )(Ca 2+ ) (CO32- )(CaCO30 )

K1 =

17

(Na + ) (CO32- )K8 =(NaCO3- )(Na + ) (SO 2-4 )K7 =(NaSO-4 )K10 =(K + ) (SO 2-4 )(KSO-4 )(H + ) (H 3SiO-4 )K11 =(H 4SiO 4 )K12 =(Na + ) (HCO3- )K9 =(NaHCO30 )(H + ) 2 (H 2SiO 2-4 )(H 4SiO 4 )(3.11)(3.12)(3.13)

where Ki are the equilibrium constants of the ith complexed species [-].

3.5. Cation Exchange and Selectivity
Partitioning between the solid and solution phases is described with the Gapon equation [Whiteand Zelazny, 1986]
ciy + ( c j )x + 1/ x

K ij =

c xj + ( ciy + )1/ y(3.14)

where y and x are the valences of species i and j, respectively, and Kij is the Gapon selectivity
coefficient [-]. The adsorbed concentration is expressed in (molckg-1 soil). It is assumed that the
cation exchange capacity cT (molckg-1 soil) is constant and independent of pH.
cT = ∑ ci(3.15)

When four cations ( Ca , Mg , Na and K ) are involved in the exchange reactions, the following
system of equations results:cT = Ca 2+ + Mg 2+ + Na + + K +)Ca ( Mg )( Na )Mg
K =Na ( Mg )
(K )
Mg
K =
K ( Mg )
K13 =

14

15

Mg 2+ ( Ca

2 + 1/ 22+

2 + 1/ 22+++

2 + 1/ 22+++

2 + 1/ 218(3.16)(3.17)

3.6. Precipitation-Dissolution Reactions
The carbonate chemistry module considers four solid phases which, if specified or approachedfrom oversaturation, must be in equilibrium with the solution: gypsum, nesquehonite,hydromagnesite, and sepiolite. Precipitation-dissolution of calcite can be optionally treatedassuming either equilibrium or by means of rate equations. In the latter case, the equation
corresponding to calcite equilibrium presented in this section is omitted from the equilibrium
system, and the rate of calcite precipitation-dissolution is calculated from the rate equation as
described later. Dissolution of dolomite, which will also be discussed later, is always considered
as a kinetic process and never included in an equilibrium system since ordered dolomite almost
never precipitates in soils. We refer to Suarez and Šimůnek [1997] for a detailed discussion on
how to select and consider these solids. The precipitation or dissolution of gypsum, calcite (ifconsidered in the equilibrium system), nesquehonite, hydromagnesite, and sepiolite in thepresence of CO2 are described by
CaSO 4 • 2H 2 O ⇔ Ca 2+ + SO 42− + 2H 2 O(3.18)

CaCO3 + CO 2 (g) + H 2 O ⇔ Ca 2+ + 2HCO3−(3.19)

MgCO3 • 3H 2 O ⇔ Ca 2+ + 2HCO3− + 2H 2 O(3.20)

Mg 5 (CO3 ) 4 (OH) 2 • 4H 2 O + 6 CO 2 (g) ⇔ 5Mg 2+ + 10HCO3−(3.21)

Mg 2Si3O7.5 (OH) • 3H 2 O + 4.5 H 2 O + 4CO 2 (g) ⇔ 2Mg 2+ + 3H 4SiO 4 + 4HCO3−(3.22)

with the corresponding solubility products KSP [-] given byG
K SP= (Ca 2+ ) (SO 24− ) (H 2 O) 2

(3.23)C
K SP= (Ca 2+ ) (CO32− )(3.24)N
K SP= (Mg 2+ ) (SO 24− ) (H 2 O)3(3.25)H
K SP= (Mg 2+ )5 (CO32− ) (OH - ) 2 (H 2 O) 4

(3.26)S
K SP=(Mg 2+ ) 2 (H 4SiO 4 )3 (OH - ) 4(H 2 O) 4.5

19(3.27)

where the indexes G, C, N, H, and S refer to gypsum, calcite, nesquehonite, hydromagnesite andsepiolite, respectively.
Substituting (3.5) through (3.8) into (3.23) through (3.27) gives the solubility products for the
carbonate solids expressed in terms of bicarbonate, which is almost always the major carbonate
ion for conditions (6<pH<10.5) in which the carbonate chemistry module is assumed to be
applicable:
C(Ca 2+ ) (HCO3− ) 2 = K SP

K CO2 K a1 PCO2 (H 2 O)
Ka2
K CO2 K a1 PCO2
N

(Mg 2+ ) (HCO3− ) 2 = K SP

H
(Mg 2+ )5 (HCO3− )10 = K SP

2+ 2

− 4
3

(Mg ) (HCO ) = K

S
SP

Ka2

(H 2 O) 26
6
K CO
K a61 PCO2
2

K a42 K w2

4
4
K CO
K a41 PCO
(H 2 O) 4.5
2
2

K w4 (H 4SiO 4 )3

(3.28)(3.29)(3.30)(3.31)

Expressing the solubility products in this way significantly decreases the number of numericaliterations necessary to reach equilibrium as compared to when equations (3.23) through (3.27)
are used.

3.7. Kinetic Model for Calcite Precipitation-Dissolution
The reaction rates of calcite precipitation-dissolution in the absence of inhibitors such as "foreignions" and dissolved organic matter, RC (mmol cm-2s-1), were calculated with the rate equation of
Plummer et al. [1978]
R C =k1 (H + ) + k2 (H 2 CO*3 ) + k3 (H 2 O) − k4

Ka2
(Ca 2+ ) (HCO3− )C
K SP(3.32)

where
1
k4 =
k1 + +  k2 (H 2 CO*3 ) + k3 (H 2 O) (H s )(3.33)

and where k1, k2, and k3 are temperature-dependent first-order rate constants representing theforward reactions (mmol cm-2s-1), and k4 is a function dependent on both temperature and CO2
concentration representing the back reactions (mmol cm-2s-1). The precipitation-dissolution rate
RC is expressed in mmol of calcite per cm2 of surface area per second. The term (HS+) is the H+
activity at the calcite surface. Its value is assumed to be (H+) at calcite saturation, where20activities of H2CO3* and H2O at the calcite surface are equal to their bulk fluid values [Plummer
et al., 1978]. The temperature dependency of the constants k1, k2, k3 is expressed aslog k = a1 +a2T

(3.34)

where values of the empirical constants a1 and a2 are given by Plummer et al. [1978]. For
conditions where pH>8 and pCO2<1000 Pa, an alternative expression for the precipitation rate isused, which is considered more accurate for those conditions [Inskeep and Bloom, 1985]
C
R C = - 11.82 [(Ca 2+ ) (CO32− ) - K SP]

(3.35)

with an apparent Arrhenius activation energy of 48.1 kJ mol-1 for the precipitation rate constant[Inskeep and Bloom, 1985].
The precipitation or dissolution rate of calcite is reduced by the presence of various inhibitors.
Suarez and Šimůnek [1997] developed the following function for the reduction of the
precipitation-dissolution rates due to surface poisoning by dissolved organic carbon, based on
experimental data from Inskeep and Bloom [1986]
r = exp(- b1 x - b2 x 2 - b3 x 0.5 )(3.36)

where r is the reduction constant [-], x is dissolved organic carbon (μmol l-1) and b1, b2, and b3
are regression coefficients (0.005104, 0.000426, 0.069111, respectively).

3.8. Kinetic Model of Dolomite Dissolution
The reaction rates of dolomite dissolution, RD (mmol cm-2s-1), were calculated with the rateequation of Busenberg and Plummer [1982]
R D =k1 (H + )0.5 + k2 (H 2 CO*3 )0.5 + k3 (H 2 O)0.5 − k4 (HCO3− )(3.37)

where the temperature-dependent first-order rate constants k1, k2, k3 (mmol cm-2s-1), representing
the forward reactions, and k4 (mmol cm-2s-1), representing the back reaction, are given by (3.34)
with empirical constants a1 and a2 given by Busenberg and Plummer [1982]. The dissolution rate
RD is again expressed in mmol of dolomite per cm2 of surface area per second. These rate
constants are used for ion activity products IAPD<10-19. For values below 10-19, the rate isextremely low and assumed to be zero [Busenberg and Plummer, 1982] in the absence ofadditional data.

21

3.9. Silica Concentration in Soil Solution
Relatively little information exists about Si concentrations in soil water. The use of equilibrium
calculations of silica solubility from the stable mineral (quartz) results in the unrealistic
prediction that solution concentrations are independent of pH up to pH 8, above which thesolubility will increase due to the dissociation of silicic acid. Si concentrations in soils areusually not fixed by quartz solubility but rather by dissolution (and possibly precipitation) ofaluminosilicates, including poorly crystallized phases and Si adsorption-desorption onto oxides
and aluminosilicates. As a result of these reactions, Si concentrations in soil solutions follow a
U-shaped curve with pH, similar to Al oxide solubility, with a Si minimum around pH 7.5[Suarez, 1977]. Suarez [1977] developed a simple relation between silica content in the soil
solution and the soil pH:

SiO 4T =+d1 d 2 p H + d3 p H

2(3.38)

in which the empirical constants d1, d2, and d3 are equal to 6340, 1430, and 81.9, respectively,
and where SiO4T is the sum of all silica species expressed in mol l-1. We utilize this expression
and the dissociation expressions for K11 and K12 (eq. (3.13)) only to obtain estimates of H4SiO4from total SiO4. As a result, sepiolite reactions are not expressed in terms of H3SiO4- andH2SiO42-, which are not included in the charge balance expressions. Only the species H4SiO4 isused in the carbonate chemistry module.

3.10. Activity Coefficients
3.10.1. Extended Debye-Hückel Expression
Activity coefficients in the dilute to moderately saline solution range are calculated using theextended version of the Debye-Hückel equation [Truesdell and Jones, 1974]:
ln γ = -Az 2 I+ bI
1 + Ba I(3.39)

where A (kg0.5mol-0.5) and B (kg0.5cm-1mol-0.5) are constants depending only upon the dielectric
constant, density, and temperature, z is the ionic charge in protonic units, a (cm) and b (kg mol-1)
are two adjustable parameters, and I is the ionic strength (mol kg-1):
M

I = 0.5∑ zi2 ci
i=122(3.40)

where M is the number of species in the solution mixture. The adjustable parameters a and b forindividual species are given by Truesdell and Jones [1974]. The activities of neutral species arecalculated asln γ = a′ I

(3.41)

where a' is an empirical parameter.
If the extended Debye-Hückel theory is used to calculate activity coefficients, the activity ofwater is then calculated in the same way as in the WATEQ program [Truesdell and Jones, 1974]
using the approximate relationM(H 2 O) = 1 − 0.017 ∑mi(3.42)

i =1

3.10.2. Pitzer Expressions
At high ionic strength, activity coefficients are no longer universal functions of ionic strength,
but depend on the relative concentration of the various ions present in the solution [Felmy andWeare, 1986]. The activity coefficients can then be expressed in terms of a virial-type expansionof the form [Pitzer, 1979]ln γ i = ln γ iDH + ∑Bij ( I )m j + ∑∑Cijk m j mk + ...
jj(3.43)kwhere γiDH is a modified Debye-Hückel activity coefficient which is a universal function of ionicstrength, and Bij and Cijk are specific coefficients for each interaction. The subroutines for the
calculation of the Pitzer activity coefficients were adopted from the GMIN code [Felmy, 1990].
This model is considered accurate also for solutions with very high ionic strength (up to 20 molkg-1), and can be used down to infinite dilution.
If the Pitzer theory is used, then the activity of water is obtained from the expression [Felmy andWeare, 1986]ln (H 2 O) = -W M  ∑mi  φ1000  i =1 

(3.44)

where W is the molecular weight of water and φ the osmotic coefficient (defined in Section 3.12).
We refer to Felmy and Weare [1986] for a more detailed discussion.

23

3.11. Temperature Effects
Most thermodynamic equilibrium constants depend on both the temperature and pressure of thesystem. The temperature dependence of the thermodynamic equilibrium constants is often
expressed as a power function of the absolute temperature:log K = a1 +aa2
+ a3T + a4 log T + 52
TT

(3.45)

where T is absolute temperature [K], and a1 through a5 are empirical constants. The pressure
dependence can be neglected for most soils. The empirical constants for the temperaturedependent thermodynamic constants used in the calculations are listed in Šimůnek et al. [1996].
The temperature dependence of the equilibrium constants for which the constants of equation(3.45) do not exist is expressed with the enthalpy of reaction and the Van 't Hoff expression[Truesdell and Jones, 1974].

3.12. Osmotic Coefficient and Osmotic Pressure Head
We use the semiempirical equation of Pitzer [1973] and co-workers to calculate the osmotic
coefficient φ. The osmotic pressure of electrolyte solutions, Pφ (Pa), is related to the osmotic
coefficient φ and molality as follows [Stokes, 1979]

Pφ = RTM s ν mφ
Vs m0

(3.46)

where Vs is the partial molar volume of the solvent (cm3mol-1), m0 is unit molality (1 mol kg-1),
and Ms is the molar weight (mol-1). The osmotic pressure head, hφ [L], is related to the osmoticpressure by

hφ =

Pφ

ρg

(3.47)

where ρ is the density of water [ML-3] and g is the gravitational constant [L2T-1].
3.13. Effect of Solution Composition on Hydraulic Conductivity
The accumulation of monovalent cations, such as sodium and potassium, often leads to clay
dispersion, swelling, flocculation, and overall poor soil physico-mechanical properties. These
processes have an adverse effect on the soil hydraulic properties, including hydraulicconductivity, infiltration rates, and soil water retention resulting from swelling and clay24

dispersion. These negative effects are usually explained based on the diffuse double layer theory.
A consequence of the more diffuse double layer in the presence of monovalent ions as compared
to divalent ions is the greater repulsion force or swelling pressure between neighboring clay
platelets. These negative effects become more pronounced with decreasing salt concentration andvalence of the adsorbed ions [Shainberg and Levy, 1992]. In addition, Suarez et al. [1984]
determined that elevated levels of pH also had an adverse effect on the saturated hydraulicconductivity.
The effect of solution chemistry on the hydraulic conductivity in the major ion chemistry module
is calculated as follows
K ( h, pH , SAR, C0 ) = r ( pH , SAR, C0 ) K ( h )(3.48)

where SAR is the sodium adsorption ratio, C0 is the total salt concentration of the ambient
solution in mmolcl-1, and r is a scaling factor that represents the effect of the solution
composition on the final hydraulic conductivity [-], and which is related to pH, SAR, and salinity.
The hydraulic conductivity without the scaling factor r can be assumed to be the optimal value
under favorable chemical conditions in terms of optimal pH, SAR, and salinity. Although soil
specific, the effects of solution chemistry are too important to ignore. We include reduction
functions calculated for some illitic soils of California based on the experimental work of
McNeal [1968] and Suarez et al. [1984]. The overall scaling factor r in equation (3.48) for thispurpose is divided into two partsr ( pH , SAR, C0 ) = r1 ( SAR, C0 ) r2 ( pH )(3.49)

where the first part, r1 [-], reflects the effect of the exchangeable sodium percentage and dilutionof the solution on the hydraulic conductivity, while the second part, r2 [-], represents the effect of
the soil solution pH. The first term is based on a simple clay-swelling model, which treats mixedions clays as a simple mixture of homoionic sodium and calcium clay. Clay swelling is thenrelated to a decrease in soil hydraulic conductivity [McNeal, 1974]. The r1 term was defined by
McNeal [1968] as
r1 = 1 -cx n
1 + cx n(3.50)

where c and n are empirical parameters, and x is a swelling factor. The interlayer swelling of soil
montmorillonite, x, is defined in the following way
x = f mont 3.6 •10−4 ESP*d *(3.51)

where fmont is the weight fraction of montmorillonite in the soil, d* is the adjusted interlayer
spacing [L] and ESP* is the adjusted exchangeable sodium percentage. For most soils, one can25use the assumption that fmont = 0.1 [McNeal, 1968]. The adjusted exchangeable sodium
percentage is calculated as=
ESP* max[0, ESP-(1.24 + 11.63 log C0 )](3.52)

where C0 is total salt concentration of the ambient solution in mmolcl-1, and ESP is defined as

ESP = 100

Na
CEC

(3.53)

where CEC is the soil cation exchange capacity (mmolckg-1) and Na the exchangeable sodium
concentration (mmolckg-1). The adjusted interlayer spacing, d*, is given by

for C0 > 300 mmolc liter −1

*

d =0
*

−1/ 2

d = 356.4 C0+1.2

for C0 < 300 mmolc liter −1(3.54)

McNeal [1968] reported that the values of the empirical factor n in equation (3.50) depends
primarily on the soil ESP and that as a first approximation n values may be estimated using

n =1
n=2

for ESP < 25
for 25 ≤ ESP ≤ 50

n=3

for ESP > 50

(3.55)

Only the values of the empirical factor c vary from soil to soil. HYDRUS uses values reported by
McNeal [1968]:

c = 35c = 932

for ESP < 25
for 25 ≤ ESP≤ 50

c = 25 000

for ESP > 50

(3.56)The reduction factor, r2, for the effect of pH on hydraulic conductivity was calculated fromexperimental data of Suarez et al. [1984] after first correcting for the adverse effects of low
salinity and high exchangeable sodium using the r1 values. The following equation was used
r2 = 1
r2 = 3.46 - 0.36 p H
r2 = 0.1

for p H< 6.83
for 6.83 ≤ p H≤ 9.3
for p H > 9.3

26(3.57)

Note that although the models for reductions in the hydraulic conductivity due to changes in thesolution composition were derived from data on the saturated hydraulic conductivity, the same
reduction factors were used for the entire range of the pressure heads. The assumption that the rvalues for saturated conditions can be applied to the entire range of pressure heads has not beenclosely examined thus far in the literature.

27

28

4. Preprocessing
4.1. Main Processes
The UNSATCHEM module is activated in the Main Processes window (Fig. 1) by selecting the
Solute Transport check box and the Major Ion Chemistry (Unsatchem) radio button.

Figure 4.1. The Main Processes dialog window for the UNSATCHEM module.

4.2. Solute Transport - General Information
When the UNSATCHEM module is used, basic information needed for defining the solutetransport problem is entered in the Solute Transport dialog window displayed in Figure 4.2. In
this window, users specify again the Space and Time Weighting Schemes and additional Solute
Information such as mass units. Note that all concentrations in the UNSATCHEM module are
either in meq/L (in the liquid phase) or meq/kg (in the solid phase) (meq=mmolc). The pulse
duration is not used in UNSATCHEM, and the number of solutes is fixed to 8 (i.e., the numberof considered major ions: Ca2+, Mg2+, Na+, K+, Alkalinity, SO42-, Cl-, and an independent tracer).
The number of solution, adsorbed, and precipitated concentration combinations is specified whensimulating the transport of major ions in the UNSATCHEM module. This value represents the
maximum number of (solution, adsorbed, and precipitated) concentration combinations, which
can be used to specify the initial and boundary (e.g., solution composition of rainwater, waterused for irrigation, groundwater, etc.) conditions.

29

Figure 4.2. The Solute Transport dialog window for the UNSATCHEM module.

4.3. Solution Compositions for the UNSATCHEM Module
The set of solution, adsorbed, and precipitated concentration combinations for the
UNSATCHEM module is specified in the Solution Compositions dialog window (Fig. 4.3).
These (solution, adsorbed, and precipitated) concentration combinations can be used to specify
the initial and boundary conditions. The number of solution, adsorbed, and precipitatedconcentration combinations is specified in the General Solute Transport Information window(Fig. 4.2). Solution Concentrations need to be specified for all major ions: Ca2+, Mg2+, Na+, K+,
Alkalinity, SO42-, Cl-, and an independent tracer; Adsorbed Concentrations for all cations: Ca2+,
Mg2+, Na+, and K+; and Precipitated Concentrations for all solids that UNSATCHEM canconsider: calcite, dolomite, gypsum, nesquohonite, hydromagnesite, and sepiolite. Solution
Concentrations need to be specified in meq/L (L = liter), and Adsorbed Concentrations andPrecipitated Concentrations in meq/kg (meq=mmolc).

30

Figure 4.3. The Solution Compositions dialog window for the UNSATCHEM module.

4.4. Solute Reaction Parameters
When the UNSATCHEM module is used, the Solute Reaction (and transport) Parameters are
specified in the Solute Reaction Parameters dialog window displayed in Figure 4.4. The
following Soil Specific Parameters are specified for each soil material:
Bulk.d.
Dw
Disper.L.
Disper.T.
CEC
Calc.SA
Dol.SA
DOC
K[Ca/Mg]
K[Ca/Na]K[Ca/K]

Bulk density, ρ [g/cm3]
Molecular diffusion coefficient in free water, Dw [L2T-1]
Longitudinal dispersivity, DL [L]Transverse dispersivity, DT [L]
Cation exchange capacity, CEC [meq/kg]
Calcite surface area [m2/dm3]
Dolomite surface area [m2/dm3]
Dissolved organic carbon [µmol/dm3]
Gapon constant for exchange of calcium and magnesium
Gapon constant for exchange of calcium and sodium
Gapon constant for exchange of calcium and potassium

31

Figure 4.4. The Solute Reaction Parameters dialog window for the UNSATCHEM module.

4.5. Chemical Parameters for the UNSATCHEM module.
The following chemical parameters and selections for the UNSATCHEM module are specified
in the Chemical Parameters dialog window (Fig. 4.5):- Whether the kinetic or equilibrium model for the precipitation and dissolution of calcite
and dissolution of dolomite is to be used (the Kinetic Precipitation/Dissolution checkbox).
- Whether the silica content in the solution is to be calculated based on the solution pH orwhether the effect of pH is to be neglected (the Silica in Solution (pH Dependency) checkbox).
- The Critical Ionic Strength, i.e., the ionic strength, below which the extended DebyeHückel equation is used to calculate ion activity coefficients. Pitzer virial-type equations
are used above this value.
- The Maximum Number of Iterations allowed during any time step between the solutetransport and chemical modules. When the maximum number of iterations is reached, thenthe code proceeds to the new time level. The recommended value (if the iterative approachis to be used; from our experience) is 5. Set equal to one if no iteration (we recommend thisnon-iterative approach) is required (this, in general, leads to significantly lower
computational time without significantly altering the results in most cases).
- Whether the hydraulic conductivity is modified depending on the solution chemistry using
the McNeal [1968] semiempirical approach (the Conductivity Reduction due toChemistry check box).

32

Figure 4.5. The Chemical Parameters dialog window for the UNSATCHEM module.

4.6. Initial Conditions
The initial conditions for the UNSATCHEM module are defined in terms of the Solution
Composition numbers (integer), Exchange Species numbers, Solid Species numbers, and CO2
Concentrations (real). The composition numbers refer to different Solution Compositions,Exchange Species, and Solid Species defined in Figure 4.3.

Figure 4.5. The "Initial Conditions" part of the Data Tab of the Navigator Bar for the
UNSATCHEM module.

33

34

5. Post-processing
5.1. Results – Graphical Display
Multiple variables can be displayed in the View window (Fig. 5.1). A comprehensive list ofvariables that can be displayed in the View window is summarized in Table 5.1.

Table 5.1. Additional UNSATCHEM variables displayed in the View Window of the Results tab(Results - Graphical Display).
--------------------------------------------------------------------------------------------------------------------Variable
Description--------------------------------------------------------------------------------------------------------------------UNSATCHEM Variables Major ions (Ca, Mg, Na, K, HCO3, SO4, Cl [in mmolcl-1], tracer [-],
sorbed Ca, sorbed Mg, sorbed Na, sorbed K, calcite, gypsum,
dolomite, nesquohonite, hydromagnesite, sepiolite (in mmolckg-1),
pH [-], EC [dS/m], and SAR.
---------------------------------------------------------------------------------------------------------------------

Figure 5.1. The "Results - Graphical Display" part of the Data Tab of the Navigator Bar for the
UNSATCHEM module.

355.2. Results – Other Information
The Chemical Mass Balance information (Fig. 5.2) (on the Data Tab of the Navigator Bar or at
the Results Menu) is provided for the UNSATCHEM module. This file gives mass balances indifferent phases (liquid, sorbed, solid) for all major ions (i.e., Ca, Mg, Na, K, HCO3, SO4, and Cl).

Figure 5.2. The Chemical Mass Balance Information dialog window for the UNSATCHEM
module.

36

6. Input and Output Files
The UNSATCHEM module requires several additional input files (described in Section 6.1
below) providing various coefficients for calculating ion activities using the Pitzer's expressions
and has a slightly different format of the Selector.in file (described in Section 6.2 below) in thesection related to solute transport.
6.1. Pitzer's Coefficients
The major ion chemistry module requires four additional input files containing input data for the
Pitzer equations. These four input files are provided together with the program and should not bechanged by the user. The four input files, which were adopted from Felmy [1990], are placed in
the ChemData directory under the executable program.
COMP.DAT

contains the species ID numbers, species names, and species charge.

BINARYP.DAT

contains the ID number of each species in each binary interaction
considered (e.g., CaHCO3+) and the Pitzer ion interaction parameters β(0),β(1), β(2), and Cφ for binary systems.

TERNARYP.DAT

contains the Pitzer ion-interaction parameters for the common ion ternary
systems, θ, and ψ. The first two columns include the cation-cation oranion-anion ID numbers associated with the ion-interaction parameter, θ,in column three. Subsequent columns include the anion or cation ID
number and the triple ion-interaction parameter, ψ, associated with that
triple ion interaction.

LAMBDA.DAT

contains the Pitzer ion-interaction parameters for the neutral species, λand δ. The first column of this file contains the ID number for the neutral
species, and the second column the ID number for the cation or anion
involved in the neutral-cation or neutral-anion interaction parameterized
by the Pitzer λ parameter included in the third column. Subsequentcolumns are for higher-order neutral interactions.

6.2. The Selector.in Input FileHYDRUS (2D/3D) with UNSATCHEM module has a different structure of the part of theSelector.in file that is related to solute transport. This is described in detail in Table 6.1.

37

Table. 6.1. Major ion chemistry information in the Selector.in input file.
RecordType

Variable

Description

1,2

-

-

Comment lines.

3

Real

Epsi

Temporal weighing coefficient.
=0.0 for an explicit scheme,=0.5 for a Crank-Nicholson implicit scheme.
=1.0 for a fully implicit scheme.

3

LogicallUpW.true. if upstream weighing formulation is to be used (see Section 6.4.2 of the
HYDRUS technical manual).
.false. if the original Galerkin formulation is to be used.

3

LogicallArtD.true. if artificial dispersion is to be added in order to fulfill the PeCr stability
criterion (see Section 6.4.6 of the HYDRUS technical manual).
.false. otherwise.

3

Real

PeCr

Stability criteria (see Section 6.4.6 of the HYDRUS technical manual). Set equal
to zero when lUpW is equal to .true..

4

-

-

Comment line.

5

LogicallRate

Set this logical variable equal to .true. if kinetic precipitation-dissolution of
calcite and kinetic dissolution of dolomite is to be considered.
Set this logical variable equal to .false. if only equilibrium reactions are to beconsidered.

5

LogicallSilicSet this logical variable equal to .true. if the silica content of the soil solution is
to be calculated based on the solution pH.
Set this logical variable equal to .false. if the silica content of the soil solution is
not considered.

5

Real

UCrit

Value of ionic strength below which the extended Debye-Hückel equation isused to calculate ion activity coefficients. Pitzer's virial-type equations are used
above this value. It is suggested that either one or the other be used.

5

Integer

MaxCh

Maximum number of iterations allowed during any time step between the solutetransport and chemical modules. When the maximum number of iterations isreached for MaxCh>5 then the time step is divided by three and the particulartime level is restarted again. When the maximum number of iterations is reached
for MaxCh≤5 then the code proceeds to the new time level. Recommended value(from experience) is 5. Set equal to one if no iteration is required.

5

LogicallkRed

Set this logical variable equal to .true. if hydraulic conductivity is to be adjusted
depending on the solution chemistry (see Section 3.13).
Set this logical variable equal to .false. if hydraulic conductivity is not to be
adjusted depending on the solution chemistry (see Section 3.13)..

6

-

-

Comment line.

7
7
7
7
7
7

Real
Real
Real
Real
Real
Real

ChPar(1,M) Bulk density of material M [g/cm3].
ChPar(2,M) Molecular diffusion coefficient in free water, Dw [L2T-1].
ChPar(3,M) Longitudinal dispersivity for material type M, DL [L].
ChPar(4,M) Transverse dispersivity for material type M, DT [L].
ChPar(5,M) Cation exchange capacity for material type M, cT (mmolc kg-1 of soil).
ChPar(6,M) Calcite surface area AC (m2l-1 of soil matrix). Set this variable equal to zero
when lRate=.false..

38

Table 6.1. (continued)

RecordType

Variable

Description

7

Real

ChPar(7,M)

7

Real

ChPar(8,M)

7
7
7

Real
Real
Real

ChPar(9,M)
ChPar(10,M)
ChPar(11,M)

Dolomite surface area AD (m2l-1 of soil matrix). Set this variable equal to zero
when lRate=.false. or dolomite is not present in the soil profile.
Dissolved organic carbon (μmol l-1). This variable is used to calculate thereduction in the precipitation-dissolution rates of calcite.
Gapon' selectivity constant between calcium and magnesium, K13 [-].
Gapon' selectivity constant between calcium and sodium, K14 [-].
Gapon' selectivity constant between calcium and potassium, K15 [-].
Record 7 information is provided for each material M (from 1 to NMat).

8

-

-

Comment line.

9

IntegernSolConc

Number of different solutions used in a particular application.

9

IntegernAdsConc

Number of surface species combinations used in a particular application.

9

IntegernPrecConc

Number of mineral phase combinations used in a particular application.

10

-

-

Comment line.

11

Real

ConcTab(1,1)

11

Real

ConcTab(1,2)

11

Real

ConcTab(1,3)

11

Real

ConcTab(1,4)

11

Real

ConcTab(1,5)

11

Real

ConcTab(1,6)

11

Real

ConcTab(1,7)

11

Real

ConcTab(1,8)

Analytical concentration of calcium for the first solution concentration
combination, CaT (mmolcl-1 of solution).
Analytical concentration of magnesium for the first solution concentration
combination, MgT (mmolcl-1 of solution).
Analytical concentration of sodium for the first solution concentration
combination, NaT (mmolcl-1 of solution).
Analytical concentration of potassium for the first solution concentration
combination, KT (mmolcl-1 of solution).
Analytical concentration of alkalinity for the first solution concentration
combination, Alk (mmolcl-1 of solution).
Analytical concentration of sulfate for the first solution concentration
combination, SO4T (mmolcl-1 of solution).
Analytical concentration of chloride for the first solution concentration
combination, ClT (mmolcl-1 of solution).
Analytical concentration of hypothetical tracer for the first solution
concentration combination [-].
In general, one record as described above is required for each solution
concentration combination, starting with the first solution concentration
combination and continuing in sequence up to the nSolConcth combination.

12

-

-

Comment line.

13

Real

XConcTab(1,1) Adsorbed (surface species) calcium concentration for the first surface
species combination, a (mmolckg-1 of soil matrix).

13

Real

XConcTab(1,2) Adsorbed magnesium concentration for
combination, g (mmolckg-1 of soil matrix).

13

Real

XConcTab(1,3) Adsorbed sodium concentration for the first surface species combination, a(mmolckg-1 of soil matrix).

13

Real

XConcTab(1,4) Adsorbed potassium concentration for the surface species combination,(mmolckg-1 of soil matrix).

39

the

first

surfacespeciesTable 6.1. (continued)

RecordType

Variable

Description
In general, one record as described above is required for each surface species
combination, starting with the first surface species combination andcontinuing in sequence up to the nAdsConcth combination.

14

-

-

Comment line.

15

Real

SConcTab(1,1) Solid-phase calcite concentration for the first mineral phase combination
expressed in mmolc of Ca per kg of the soil matrix, CaCO3 (divide by 2*10-3
to obtain moles of calcite per kg of soil matrix).

17

Real

SConcTab(1,2) Solid-phase gypsum concentration for the first mineral phase combination
expressed in mmolc of Ca per kg of the soil matrix, CaSO4 (divide by 2*10-3
to obtain moles of gypsum per kg of soil matrix).

17

Real

SConcTab(1,3) Solid-phase dolomite concentration for the first mineral phase combination
expressed in mmolc of Ca per kg of the soil matrix, CaMg(CO3)2 (divide by
2*10-3 to obtain moles of dolomite per kg of soil matrix).

17

Real

SConcTab(1,4) Solid-phase hydromagnesite concentration for the first mineral phasecombination expressed in mmolc of Mg per kg of the soil matrix,
Mg5(CO3)4(OH)2· 4H2O (divide by 10-4 to obtain moles of hydromagnesite
per kg of soil matrix).

17

Real

SConcTab(1,5) Solid-phase nesquehonite concentration for the first mineral phasecombination expressed in mmolc of Mg per kg of the soil matrix,
MgCO3·3H2O (divide by 2*10-3 to obtain moles of nesquehonite per kg ofsoil matrix).

17

Real

SConcTab(1,6) Solid-phase sepiolite concentration for the first mineral phase combination
expressed in mmolc of Mg per kg of the soil matrix, Mg2Si3O7.5(OH)·3H2O(divide by 4*10-3 to obtain moles of nesquehonite per kg of soil matrix).
In general, one record as described above is required for each mineral phasecombination, starting with the first mineral phase combination andcontinuing in sequence up to the nPrecConcth combination.

40

6.3. UNSATCHEM Output Files
The UNSATCHEM module prints at prescribed print times, in addition to output printed by thestandard HYDRUS module, the following files:
CONCi.OUT

Nodal values of the aqueous concentrations of calcium, magnesium,
sodium, potassium, alkalinity, sulfate, chloride (mmolcl-1), andhypothetical tracer [-] (where i=1 through 8)

SORBi.OUT

Nodal values of the mineral phase concentrations for calcite, gypsum,
dolomite, hydromagnesite, nesquehonite, and sepiolite, and nodal valuesof the adsorbed concentrations for calcium, magnesium, sodium, andpotassium (mmolckg-1) (where i=1 through 10). Adsorbed concentrations
are printed only when cation exchange is considered.

EQUIL.OUT

This file contains the chemical information such as activities of calcium,
bicarbonate, and water, alkalinity, pH, SAR, electric conductivity of thesolution, ionic strength, osmotic coefficient, osmotic pressure head, andion activity products for calcite, gypsum, and dolomite (see Table 6.2).

CHEMBAL.OUT

This file contains the information about the total amount of particularspecies (e.g. Ca, Mg, SO4,..) in solution, mineral phase, and surface
species form in the entire flow region, as well as the cumulative
boundary fluxes and absolute mass error in particular species. All massbalances and cumulative solute fluxes are expressed in units of meq/m or
meq for two- or three-dimensional applications, respectively.

PH.OUT

Nodal values of pH.

SAR.OUT

Nodal values of Sodium Adsorption Ratio (SAR).

EC.OUT

Nodal values of Electric Conductivity (EC) (dSm-1).

41

Table 6.2. EQUIL.OUT - chemical information.

Node

Number of nodal point n.

Widthx-coordinate of node n [L].

Depthy-coordinate of node n [L].

aCa

Activity of Ca2+ [-].

aHCO3

Activity of HCO3- [-].

aH2O

Activity of water [-].

Alk

Alkalinity (mmolckg-1).

pH

Negative logarithm of hydrogen activity, -log(H), [-].

SAR

Sodium adsorption ratio, defined as [Na/(Ca+Mg)0.5] (mmol0.5l-0.5).

ECElectric conductivity of the soil solution (dSm-1) (calculated using method 3 of
McNeal et al., 1970).

U

Ionic strength (mol kg-1).

pIAP(c)Negative logarithm of the ion activity product for calcite, -log[(Ca2+)(CO32-)], [-].

pIAP(g)Negative logarithm of the ion activity product for gypsum, -log[(Ca2+)(SO42-)(H2O)2], [-].

pIAP(d)Negative logarithm of the ion activity product for dolomite, -log[(Ca2+) (Mg2+)(CO32-)2], [-].

phi

Osmotic coefficient [-].

hphi

Osmotic pressure head [m].

42

7. Example Problems
The module UNSATCHEM is developed from the variably saturated solute transport model, andthus the water flow and solute transport parts of the model have been tested earlier. Therefore, both
example problems solved in this section concentrate on demonstrating the chemical features of themodel.
Two example problems are presented in this section. The first example simulates an infiltration
process into a relatively dry soil column. This example shows predicted differences between theequilibrium and kinetic precipitation-dissolution model. The second example simulates a furrow
irrigation problem
7.1. Example 1 - Column Infiltration
This example is based on the one used in the documentation of the code SWMS_2D [Šimůnek et al.,1992] to compare results obtained with the SWMS_2D and UNSAT2 [Davis and Neuman, 1983]
codes. It originally simulated a one-dimensional laboratory infiltration experiment discussed by
Skaggs et al. [1970]. This example was altered for the purpose of this report to include
multicomponent chemical transport.
The soil water retention and relative hydraulic conductivity functions of the sandy soil are θr=0.02,
θa=-0.02, θs=θm=0.35, θk=0.2875, α=0.0410 cm-1, n=1.964, Ks=0.000722 cm s-1, Kk=0.000695 cm s1). The sand was assumed to be at an initial pressure head of -150 cm. The soil hydraulic properties
were assumed to be homogenous and isotropic. The column was subjected to ponded infiltration (aDirichlet boundary condition) at the soil surface, resulting in one-dimensional vertical water flow.
The open bottom boundary of the soil column was simulated by implementing a no-flow boundary
condition during unsaturated flow (h<0) and a seepage face with h=0 when the bottom boundary
becomes saturated (this last condition was not reached during the simulation). The solution
composition of the water initially present in the soil profile is that of calcite supersaturated well
water from the Wellton-Mohawk Irrigation District (well # 15 [Suarez, 1977]: CaT=12.2, MgT=9.7,
NaT=37.5, KT=0.3, ClT=31.1, SO4T=22.1, alkalinity=6.5 mmolcL-1). Calcite undersaturated
Colorado river water from the Grand Valley ([Rhoades and Suarez, 1977]: CaT=2.63, MgT=1.05,
NaT=2.55, KT=0.06, ClT=1.94, SO4T=2.03, alkalinity=2.33 mmolcL-1) was used as the solution
composition of the infiltrating water. The simulation was run at a temperature of 25 oC, and the soil
CO2 concentration was assumed to be equal to the atmospheric CO2 (0.000333). The precipitationdissolution of calcite was considered either as an equilibrium or kinetic process with a calcitesurface area of 0.02 m2L-1 of soil. The bulk density of the soil was taken as 1.3 g cm-3 and molecular
diffusion as 0.02 cm2s-1. Longitudinal dispersivity was equal to zero. Activity coefficients were
calculated with the Debye-Hückel equation. Cation exchange was not considered in this example.
Figures 7.1 and 7.2 show the water content and tracer concentration profiles, respectively, at various
times. Calcium concentration and alkalinity profiles at multiple times for both equilibrium andkinetic calcite precipitation-dissolution are shown in Figures 7.3 and 7.4, respectively. Calcite
profiles at various times for both equilibrium and kinetic calcite precipitation-dissolution are shownin Figure 7.5.

43

0

T0
-10

T1
-20

T2

-30

T3

-40

T4

-50

T5
T6

-60
0 .0 0 0 .0 5 0 .1 0 0 .1 5 0 .2 0 0 .2 5 0 .3 0 0 .3 5

W ater Co n ten t [-]Figure 7.1. Water content profiles at various times for example 1.

0
-10
-20
-30
-40
-50
-60
0 .0 0 .1 0 .2 0 .3 0 .4 0 .5 0 .6 0 .7 0 .8 0 .9 1 .0

Co n cen tratio n [-]Figure 7.2. Tracer concentration profiles at various times for example 1.

44

0

0

-10

-10

-20

-20

-30

-30

T3

-40

-40

T4

-50

-50

T5

-60

-60

T0
T1

6

5

4

3

2

1

7

T2

T6
2

4

Calciu m [meq /L ]

6

8

10

12

14

Calciu m [meq /L ]Figure 7.3. Calcium concentration profiles at various times for equilibrium (left) and kinetic (right)
calcite precipitation-dissolution for example 1.
00Depth [m]

10

T0

-10

20

T1

0

-20

60
900

30

1800
40

2700
3600
5400

50

7200

T2

-30

T3

-40

T4
T5

-50

T6

60
0.8

0.85

0.9

0.95

1

1.05

1.1

-60
0

Alkalinity [meq/L]

1

2

3

4

5

6

7

Alkalin ity [meq /L ]Figure 7.4. Alkalinity profiles at various times for equilibrium (left) and kinetic (right) calcite
precipitation-dissolution for example 1.

45

0

0

-10

-10

-20

-20

-30

-30

-40

-40

-50

-50
-60
0 .0

-60
0

5

10

15

20

25

0 .1

0 .2

0 .3

0 .4

0 .5

0 .6

Calcite [meq /kg ]Calcite [meq /kg ]Figure 7.5. Calcite profiles at various times for equilibrium (left) and kinetic (right) calcite
precipitation-dissolution for example 1.

46

7. 2. Example 2 - Two-Dimensional Furrow Irrigation Problem
A furrow irrigation problem was used to simulate the two-dimensional infiltration of gypsum
saturated water into a sodic soil. The simulation of sodic soil reclamation demonstrates the cationexchange feature of UNSATCHEM. The schematic representation of the flow domain for theconsidered furrow irrigation and the finite element mesh are presented in Figure 7.6. It is assumed
that every other furrow is flooded with water and that the water level in the irrigated furrow is kept
constant at a level of 6 cm. Due to symmetry, it is necessary to carry out the simulation only for thedomain between the axis of two neighboring furrows. Free drainage is used as the bottom boundary
condition, and zero flux is considered on the rest of the boundary. The initial pressure head
condition is -200 cm, and the soil hydraulic properties are the same as those used in the second
example.
The calculation was run at a constant temperature of 25 oC and a CO2 concentration of 0.01 cm3cm-3.
Root water uptake and evaporation were neglected. The bulk density of the soil was taken as 1.4 g
cm-3 and molecular diffusion as 2 cm2day-1. Longitudinal and transverse dispersivities were equal to
2 and 0.2 cm, respectively.
The solution composition of the water initially present in the soil profile is that of the following
highly sodic water: CaT=0.2, MgT=0.0, NaT=4.8, KT=0.0, ClT=4.6, SO4T=0.0, Alk=0.4 mmolcL-1.
The cation exchange capacity is equal to 100 mmolckg-1 and is divided between exchangeablecalcium and sodium (a=5.0, a=95.0 mmolckg-1). The Gapon selectivity coefficients from Wagenet
and Hutson [1987] were used (K13=0.896, K14=1.158, K15=0.111). The solution composition of theirrigation water was almost gypsum saturated: CaT=32.6, MgT=0.0, NaT=4.8, KT=0.0, ClT=5.0,
SO4T=32.0, Alk=0.4 mmolcL-1. As a consequence of the reactions of the irrigation water with theexchanger composition, cation exchange was the dominant chemical processes in the soil profile.
Cation exchange is treated as an instantaneous process in the model.
Figure 7.7 shows the pressure head profiles for four different times. The distribution of the
hypothetical tracer is shown on Figure 7.8. Figures 7.9 and 7.10 present the exchangeable
concentrations of calcium and sodium, respectively. The exchange phase concentrations reflect thechanges in aqueous Na and Ca concentration. Note the significant lagging of the exchanger front toboth the water and tracer front. Also, due to the high nonlinearity of Na-Ca exchange, theconcentration and exchange fronts are very sharp, in contrast to the more diffuse tracer fronts. After
5 days, less than 25% of the profile has been reclaimed to exchangeable sodium percentage less than
15 (which has been the criteria defining the designation sodic soil). The selected cation exchange
capacity of 100 mmolckg-1 is relatively low. Selection of a higher exchange capacity and associatedhydraulic properties of a finer textured soil would enhance both the time required for infiltration as
well as the quantity of water needed for reclamation.

47

15

15

40

15

15

15

85

100 cm100 cmFigure 7.6. Schematic representation and finite element mesh of the flow domain for the furrow
irrigation system for example 2.

48

Figure 7.7. Pressure head (cm) profiles at times: a) 0.1, b) 0.5, c) 1, and d) 2 days for example 2.

49

Figure 7.8. Tracer concentration (-) profiles at times: a) 0.1, b) 1, c) 3, and d) 5 days for example 2.

50

Figure 7.9. Sodium concentration (mmolcL-1) profiles at times: a) 0.5, b) 1, c) 3, and d) 5 days forexample 2.

51

Figure 7.10. Exchangeable concentrations of sodium (mmolckg-1) profiles at times: a) 0.5, b) 1, c)
3, and d) 5 days for example 2.

52

References
Busenberg, E., and L. N. Plummer, The kinetics of dissolution of dolomite in CO2-H2O systems
at 1.5 to 65 C and 0 to 1 ATM PCO2, Am. J. Sci., 282, 45-78, 1982.
Davis, L. A., and S. P. Neuman, Documentation and user's guide: UNSAT2 - Variably saturated
flow model, Final Rep., WWL/TM-1791-1, Water, Waste & Land, Inc., Ft. Collins, CO, 1983.
Felmy, A. R., and J. H. Weare, The prediction of borate mineral equilibria in natural waters:
Application to Searles Lake, California, Geochim. Cosmochim. Acta, 50, 2771-2783, 1986.
Felmy, A. R., GMIN: A computerized chemical equilibrium model using a constrained
minimization of the Gibbs free energy, Pacific Northwest Lab., Richland, Washington, 1990.
Inskeep, W. P., and P. R. Bloom, An evaluation of rate equations for calcite precipitation kinetics
at pCO2 less than 0.01 atm and pH greater than 8, Geochim. Cosmochim. Acta, 49, 21652180, 1985.
Inskeep, W. P., and P. R. Bloom, Kinetics of calcite precipitation in the presence of watersoluble organic ligands, Soil Sci. Soc. Am. J., 50, 1167-1172, 1986.
Liu, Chen Wuing, and T. N. Narasimhan, Redox-controlled multiple-species reactive chemicaltransport. 2. Verification and application, Water Resour. Res., 25(5), 883-910, 1989.
McNeal, B. L, Prediction of the effect of mixed-salt solutions on soil hydraulic conductivity, Soil
Sci. Soc. Amer. Proc., 32, 190-193, 1968.
McNeal, B. L., J. D. Oster, and J. T. Hatcher, Calculation of electrical conductivity from solution
composition data as an aid to in-situ estimation of soil salinity, Soil Sci., 110, 405-414, 1970.
McNeal, B. L, Soil salts and their effects on water movement, In: Drainage for Agriculture,
Edited by J. van Schilfgaarde, Agronomy No 17, Am. Soc. Agr., Madison, Wi. 1974.
Pitzer, K. S., Thermodynamics of electrolytes I: Theoretical basis and general equations, J. Phys.
Chem., 77, 268-277, 1973.
Pitzer, K. S., Activity Coefficients in Electrolyte Solutions, Chap. 7, CRC Press, Boca Raton, Fl.,
1979.
Plummer, L. N., T. M. Wigley, and D. L. Parkhurst, The kinetics of calcite dissolution in CO2
systems at 5o to 60oC and 0.0 to 1.0 atm CO2, Am. J. Sci., 278, 179-216, 1978.
Plummer, L. N., and E. Busenberg, The solubilities of calcite, aragonite and vaterite in CO2-H2O
solutions between 0 and 90oC, and an evaluation of the aqueous model for the system
CaCO3-CO2-H2O, Geochim. Cosmochim. Acta, 46, 1011-1040, 1982.
Šejna, M., J. Šimůnek, and M. Th. van Genuchten, The HYDRUS Software Package for Simulating
Two- and Three-Dimensional Movement of Water, Heat, and Multiple Solutes in VariablySaturated Media, User Manual, Version 2.0, PC Progress, Prague, Czech Republic, pp. 280,
2011.
Shainberg, I., and G. J. Levy, Physico-chemical effects of salts upon infiltration and water
movement in soils, In: Interacting Processes in Soil Science, edited by R. J. Wagenet, P.
Baveye, B. A. Stewart, Lewis Publishers, CRC Press, Boca Raton, Florida, 1992.
Šimůnek, J., T. Vogel, and M. Th. van Genuchten, The SWMS_2D code for simulating water flow
and solute transport in two-dimensional variably saturated media. Version 1.1, Research Report
No. 126, U.S. Salinity Laboratory, USDA, ARS, Riverside, California, 169 pp., 1992.
Šimůnek, J., and D. L. Suarez, Two-dimensional transport model for variably saturated porous
media with major ion chemistry, Water Resour. Res., 30(4), 1115-1133, 1994.
Šimůnek, J., D. L. Suarez, and M. Šejna, The UNSATCHEM software package for simulating
one-dimensional variably saturated water flow, heat transport, carbon dioxide production

53and transport, and multicomponent solute transport with major ion equilibrium and kinetic
chemistry, Version 2.0, Research Report No. 141, U.S. Salinity Laboratory, USDA, ARS,
Riverside, California, 186 pp., 1996.
Šimůnek, J., and D. L. Suarez, Sodic soil reclamation using multicomponent transport modeling,
ASCE J. Irrig. and Drain. Engineering, 123(5), 367-376, 1997.
Šimůnek, J., M. Th. van Genuchten, and M. Šejna, Development and applications of the
HYDRUS and STANMOD software packages and related codes, Vadose Zone J 7(2), 587600, 2008.
Šimůnek, J., M. Th. van Genuchten, and M. Šejna, The HYDRUS Software Package forSimulating Two- and Three-Dimensional Movement of Water, Heat, and Multiple Solutes in
Variably-Saturated Media, Technical Manual, Version 2.0, PC Progress, Prague, Czech
Republic, pp. 258, 2011.
Skaggs, R. W., E. J. Monke, and L. F. Huggins, An approximate method for determining thehydraulic conductivity function of an unsaturated soil, Techn. Rep. No. 11, Water Resour. Res.
Center, Purdue University, Lafayette, IN, 1970.
Stokes, R. H., Thermodynamics of solutions, In: R. M. Pitkowicz, Activity Coefficients inElectrolyte Solutions, CRC Press, Inc., Boca Raton, Florida, 1979.
Suarez, D. L., Magnesium, carbonate, and silica interactions in soils, U.S. Salinity Laboratory
Annual Report, USDA, 120 pp., 1977.
Suarez, D. L., J. D. Rhoades, R. Lavado, and C. M. Grieve, Effect of pH on saturated hydraulicconductivity and soil dispersion, Soil Sci. Soc. Am. J., 48, 50-55, 1984.
Suarez, D. L., and J. Šimůnek, UNSATCHEM: Unsaturated water and solute transport modelwith equilibrium and kinetic chemistry, Soil Sci. Soc. Am. J., 61, 1633-1646, 1997.
Truesdell, A. H., and B. F. Jones, WATEQ, a computer program for calculating chemical
equilibria of natural waters, J. Res. U. S. Geol. Surv., 2(2), 233-248, 1974.
Wagenet R. J., and J. L. Hutson, LEACHM: Leaching Estimation And CHemistry Model, A
process-based model of water and solute movement, transformations, plant uptake andchemical reactions in the unsaturated zone, Continuum 2, Dept. of Agronomy, Cornell
University, Ithaca, New York, NY, 1987.
Westall, J. C., J. L. Zachary, and F. M. M. Morel, MINEQL: A computer program for the
calculation of chemical equilibrium composition of aqueous systems, Tech. Note 18, 91pp.,
Dep. of Civ. Eng., Mass. Inst. of Technol., Cambridge, 1976.
White, N., and L. W. Zelazny, Charge properties in soil colloids, In: Soil Physical Chemistry,edited by D. L. Sparks, CRC Press, BOCA Raten, Florida, 1986.
Yeh, G. T., and V. S. Tripathi, A model for simulating transport of reactive multispecies
components: Model development and demonstration, Water Resour. Res., 27(12), 30753094, 1991.

54