J. Hydrol. Hydromech., 53, 2005, 2, 73–90

ANALYSIS OF FIELD-SCALE SALINIZATION:
A COMPARISON OF DIFFERENT MODELING APPROACHES
P. J. VAUGHAN1), P.J. SHOUSE1), D.L. SUAREZ1), M. TH. VAN GENUCHTEN1*) and J.E. AYARS2)
1)

George E. Brown, Jr. Salinity Laboratory, USDA-ARS, 450 W. Big Springs Road, Riverside, CA 92507, USA;
mailto: rvang@ussl.ars.usda.gov
2)
Water Management Research Laboratory, USDA-ARS, 9611 S. Riverbend Ave., Parlier, CA 93648, USA.

Soil and groundwater salinization are major problems for irrigated agriculture in many arid and semiarid
areas of the world. Studies addressing such problems require accurate estimation of salt loadings from irrigated areas through the vadose zone to underlying groundwater. We studied Cl- transport in the vadose zone
at 45 locations in a field in the San Joaquin Valley, California, through a combination of soil sampling at six
depths (0–1.8 m) and numerical modeling using a coupled water flow and solute transport code (Unsatchem). Our purpose was to assess water and salt loadings from the heterogeneous field to groundwaterover a two-year period, and to test applicability of the code to the data. Soil sampling in November, 1995,
defined the initial water content and the Cl- concentration, and the soil hydraulic properties. Four more
sampling periods, ending in November 1997, provided data for evaluating model performance. Cl- distributions in 1997 exhibited a variety of shapes including monotonically increasing or decreasing distributions
versus depth, and profiles with maxima, minima or and sigmoidal shapes. The standard modeling approach,based on the Richards equation and the convection-dispersion equation, predicted more Cl- leaching than
was observed in the field. Somewhat improved predictions were obtained when the potential transpiration
rate was increased by a factor of 1.5. Better leaching predictions were also obtained when the model included separate mobile and immobile water fractions, mostly by improving the profile shapes. Our study
shows the importance of accurate descriptions of the lower boundary conditions, spatial variability in thewater infiltration rate, and estimation of soil surface evaporation and transpiration rates.
KEY WORDS: Solute Transport, Field Simulation, Chloride, Unsatchem.
P. J. Vaughan, P.J. Shouse, D.L. Suarez, M. Th. van Genuchten, J.E. Ayars: ANALÝZA ZASOLENIAV MIERKE POĽA: POROVNANIE ROZDIELNYCH MODELOVÝCH PRÍSTUPOV. Vodohosp. Čas.,53, 2005, 2; 40 lit., 9 obr., 3 tab.
Zasoľovanie pôdy a podzemných vôd sú hlavnými problémami pôdohospodárstva v závlahových podmienkach v mnohých arídnych a semiarídnych oblastiach sveta. Štúdie, ktoré riešia podobné problémy,vyžadujú si presné určenie zaťaženia soľami zo závlah, ktoré prechádzajú nenasýtenou oblasťou pôdy do
podzemných vôd. V tejto štúdii sa zaoberáme transportom Cl- v nenasýtenej oblasti pôdy v 45 lokalitách v
San Joaquin Valley, California, využívajúc kombináciu odberu vzoriek pôdy v šiestich hĺbkach (0–1.8 m) a
numericakým modelovaním s využitím simulačného modelu kombinovaného transportu vody a rozpustených látok (solí) (Unsatchem). Cieľom je určenie priesaku vody a rozpustených látok z heterogénnehopoľa do podzemných vôd počas dvojročného obdobia a testovanie použiteľnosti uvedeného modeluvzhľadom k vstupným hodnotám. Zo vzoriek pôdy, odobratých v novembri 1995 bolo určené počiatočné
rozdelenie vlhkosti pôdy a koncentrácia Cl-, ako aj hydraulické vlastnosti pôdy. Štyri ďalšie termíny odberov, končiac v novembri 1997, poskytli údaje pre overenie modelu. Rozdelenie Cl- v roku 1997 sa vyznačovalo rozdieľnosťou tvarov, vrátane monotónne klesajúceho alebo stúpajúceho rozdelenia koncentráciív závislosti na hĺbke, ako aj profilmi s maximom, alebo aj sigmoidálneho tvaru. Štandardný modelový prístup, založený na Richardsovej rovnici a konvektívno-disperznej rovnici, predpovedal viac vyplaveného Clako bolo pozorované v poli. O niečo lepšie výsledky boli dosiahnuté, ak sa potenciálna transpirácia zvýšila
1,5-násobne. Lepšie výsledky v prognóze vyplavovania solí boli dosiahnuté, ak model obsahoval separátne
mobilnú a imobilnú vodu, predpovšetkým zlepšením tvarov vertikálnych rozdelení koncentrácií chlóru.
V tejto štúdii bolo ukázané, aké je dôležité presné určenie dolných okrajových podmienok, priestorovej
variability rýchlosti infiltrácie, ako aj určenie výparu z povrchu pôdy a transpirácie.

73

P. J. Vaughan, P.J. Shouse, D.L. Suarez, M. Th. van Genuchten and J.E. Ayars

KĽÚČOVÉ SLOVÁ: transport rozpustených látok, simulácia v poľných podmienkach, chloridy,
Unsatchem.

1. Introduction
Soil salinization is a long-standing problem for
irrigated agriculture in arid and semiarid regions(Hillel, 1998). In the San Joaquin Valley of California the problem of soil salinization is not as
acute as in many other areas of the world, but still amajor concern, especially in areas with relatively
shallow water tables (< 3 m). Also, a positive correlation has been established between soil salinity,
boron and selenium concentrations in this area
(Deverel and Millard, 1988). Considering the soil
salinization challenges and associated environmental problems due to elevated selenium concentrations at the Kesterson reservoir resulting fromagricultural drainage (Benson et al., 1991), a needexists for accurate prediction of solute loadings
from agricultural areas that are drained.
Regional studies of salinity on the west side of
the San Joaquin Valley demonstrated that the mostsalt-affected soils are fine-grained deposits locatedon the margins of alluvial fans that are downslope
from the Diablo Range (Deverel and Gallanthine,1989). The spread of irrigated agriculture in thisarea changed the dynamics of discharge and recharge to groundwater from predominantly horizontal flow to vertical flow due to evapotranspiration and infiltration in cropped fields (Deverel andGallanthine, 1989; Belitz and Phillips, 1995). Regional studies of groundwater and its salinizationneed accurate representations of water flow and saltloadings to groundwater from the overlying vadosezone.
Variably saturated water flow and solute transport in the vadose zone are commonly modeledusing the Richards equation and the convectiondispersion equation, respectively (Nielsen et al.,
1986; van Genuchten and Simunek, 2004). Other
approaches include state-space models (Wu et al.,2001), mobile-immobile flow models (van Genuchten and Wierenga, 1976) and reformulations ofsolute transport models within a transfer function
context (Jury and Roth, 1990). Much evidence exists that the standard CDE equations often do notaccurately predict solute transport in undisturbed
field soils (Jarvis and Larsson, 2001). For example,
Dust et al. (2000) observed that near-surface bromide concentrations in a sandy soil exceeded predictions obtained with the LeachP model. Mohanty
et al. (1998) similarly were unable to predict water74

and solute transport to tile drains using traditional
flow and transport equations (Simunek et al.,1996a). Solute leaching in field soils under natural
weather conditions also could not be modeled successfully using the CDE (Snow et al., 1994;
Jacques et al., 2002). However, the CDE in several
other cases adequately represented pesticide andsolute transport under field conditions (Wagenet
and Hutson, 1986; Suarez, 2001).
Considering all of these studies, it is clear thatmodeling of vadose-zone transport in field soils issubject to considerable uncertainty. Potential
sources of uncertainty are the presence of preferential or non-equilibrium flow (Nkedi-Kizza et al.,
1985; Simunek et al., 2003) and lack of detail in thespecification of upper and lower boundary conditions, plant root water uptake, rooting distribution
and spatial variability in the material properties.
When preferential flow occurs, the flowing water
bypasses certain parts of the soil profile, thereby
resulting in more rapid transport of solutes in thepreferential flow regions and slower transport in thesoil matrix (Simunek et al., 2003). Assuming thepresence of immobile water is one approach fortreating preferential flow (van Genuchten and
Wierenga, 1976).
Despite the above uncertainties numerous modeling studies are still useful in attempts to better understand and/or predict the movement of agricultural and other contaminants in the vadose zone.
One example is the evaluation of remediation
strategies for salt-affected soils (Suarez, 2001).
Prediction of pesticide movement in soil is similarly important for development of management
practices that will maximize the effectiveness of thepesticide usage, and minimize their potential forcontamination of ground or surface waters(Wagenet and Rao, 1990).
The current study is part of a larger project designed to estimate salt loadings from irrigated areasin the San Joaquin Valley through the vadose zone
to groundwater. As part of that project we conducted extensive soil sampling within a single agricultural field during a two-year period. The objective of the current study was to test alternativemodels of chloride transport in order to evaluate
their predictive capability with respect to chloride
concentrations in the 1.8-m thick vadose zone. Agricultural productivity on the west side of the San
Joaquin Valley is affected by boron and, to a lesser

Analysis of field-scale salinization: A comparison of different modeling approaches

extent, selenium. Developing field-applicable modeling of chemical transport for Cl- is hence a necessary first step for modeling transport of these elements that form reactive oxyanions.
This study differs from previous field studies inseveral ways. First, our intention was to apply a setof recommended procedures for utilizing the models when data sets are incomplete. For example, we
utilized standard measurements of soil texture andbulk density to estimate soil hydraulic properties
based on pedotransfer functions, rather than actual
measurements. Also, we developed soil surface
boundary conditions for water flow based on the
California Irrigation Management Information System (CIMIS) estimates of the potential evapotranspiration rate, CIMIS crop coefficients, measured
precipitation, and an assumed uniform distributionof irrigation water throughout the field. Second, we
wanted to conduct simulations that ran continuously for a two-year period, thereby including
simulations of both crop growth and fallow periodsin the fall and winter. Finally, we wanted the study
to apply to the entire field with all of its naturalvariability. For this purpose we conducted soil
sampling and simulations at 45 different locationsin the field.
2. Methods
2.1 Field experiment
The study was conducted in a 65 ha agricultural
field located approximately 100 km west of Fresno
in the Broadview Water District, San Joaquin Valley, California. Soils at the field site consisted of a
Cerini series, fine-loamy, mixed (calcareous), superactive, thermic Fluventic Haplocambids, and a
Tranquillity series, very fine, smectitic, thermic
Sodic Haploxerert. The Cerini clay-loam occurredin a broad band trending SW-NE through the centerof the field. The Tranquillity clay was present in
10% of the field area in the NW and SE corners(Wu et al., 2001).
Soil sampling locations were laid out on a 6 x 6square grid with a spacing of 140.4 m. In addition
to 36 grid sites, 9 more sites at other locations
within the field were sampled (Fig. 1). These sites
were identified in 1994 during a study of local
groundwater levels. Soil sampling took place during the weeks starting November 13, 1995, April
15, 1996, July 15, 1996, May 12, 1997 and November 10, 1997. Six soil samples were collectedfrom a single borehole at each location. The coresample tube in each case was 0.05 m in diameter,with samples taken at 0.3 m increments to a depthof 1.8 m. Observation wells were located on the 6 x
6 grid and installed during the growing seasons in
1996 and 1997. These wells provided depth-towater measurements and groundwater samples forchemical analyses. The groundwater chemicalcompositions of the 9 extra sites were estimatedfrom nearby wells using inverse-distance-squared
weighting.
During the experimental study (1994 through1997) the water table in the field was manipulated
as part of a subsurface irrigation experiment. This
was accomplished by controlling flow within thefield's E-W tile-drain laterals. The seven tile-drainlaterals were spaced 123 m apart and installed at a
depth of 1.5–1.8 m. Valves installed at the easternend of tile drain laterals were adjusted to control thegroundwater level. When the valves were adjusted
for higher water levels, water pressures increased in
the tile drains and at times may have reversed thenormal flow direction.
2.2 Soil chemical analysis
Soil samples were extracted at room temperature(~22 °C) with deionized water at a 1:1 ratio basedon mass. Chloride analyses of the extract were performed using coulometric titration (Cotlove et al.,1958). The concentrations were recalculated to fieldwater contents based on conservation of Cl- mass.
2.3 Hydraulic Properties
Particle size analyses were performed on 40-g
soil samples representing six depths at each of the
45 locations. Sand fractions were determined bysieving, silt and clay fractions using the hydrometer
method. Bulk densities were determined for soilcores, 0.05 m in diameter and 0.3 m long. Prior toweighing, the soil was dried at 100°C until no further weight loss was detected. Water retention parameters, including the saturated water content (θs),
the residual water content (θr), and the shape parameters α and n (van Genuchten, 1980), were estimated using the Rosetta computer program basedon a hierarchical artificial neural network approach(Schaap et al., 1998). Estimates were obtained using the bulk density and texture as determined byparticle size analyses (Fig. 2). The saturated hydraulic conductivity (Ks) was also estimated with
Rosetta. An upper limit was established at 0.8 m d-1
75

P. J. Vaughan, P.J. Shouse, D.L. Suarez, M. Th. van Genuchten and J.E. Ayars

Fig. 1. Sampling locations and kriged surface of mean Cl- concentration in the 0.0–1.8 m depth range (November, 1997). Underlined labels signify Cl- maxima in the profile.
Obr. 1. Lokality odberu vzoriek pôdy a krigované povrchy priemernej koncentrácie Cl- v hĺbkach 0,0–1,8 m (november 1997).
Značky označujú maximá Cl- v profile.

Fig. 2. Soil texture in three fractions for 45 locations and six depths.
Obr. 2. Textúra pôdy v 45 lokalitách a v šiestich hĺbkach.

in order to obtain stable behavior in the water flow
calculations. This constraint was needed for only 2of 258 estimated values at 43 locations where computations were successful, and hence was not significant. A potentially more significant source oferror was the development of a possible capillarybarrier at locations where the hydraulic conductivity was high due to the presence of sand lenses. A
76

one-dimensional model is not suited to deal withpossible lateral flow due to this layering. Hydraulicparameters (θs, θr, α, n, Ks) were specified for six
soil layers of which the top five layers were 0.3 m
thick, while the hydraulic properties of the 1.5–1.8
m layer were extended to the lower boundary of thefinite-element model at 3 m depth.

Analysis of field-scale salinization: A comparison of different modeling approaches

2.4 Solute Transport Models
Water flow and multicomponent chemical transport were simulated for a two-year period startingin November 1995 at 45 locations in the field (Fig.
1). Other processes modeled in this simulation wereCO2 flux, heat flow, root water uptake and root
growth. The CO2 flux and heat flow, while not significant for chloride transport, were needed for acompanion study of boron transport (Vaughan etal., 2004). Water flow was calculated using a finiteelement numerical solution of the Richards equation. Solute transport was modeled using either the
standard convection-dispersion equation (CDE) or
the mobile-immobile water model (MIM). For the
CDE simulations, multicomponent chemical transport was modeled using∂θ wck∂cˆ ⎤ ∂ ⎡∂c⎡ ∂c⎤
+ ρ ⎢ k + k ⎥ = ⎢θ w D k − qck ⎥∂t∂t ⎦ ∂z ⎣
∂z
⎣ ∂t
⎦

(1)

as incorporated in the Unsatchem code (Simunek etal., 1996b). In Eq. (1), θw is the volumetric watercontent [m3 m-3], ck – the concentration of the kth
species in the soil solution [mmol L-1], ck – thetotal surface concentration of the kth species[mmol kg-1], cˆk – the total concentration of the kth
species in solid phases [mmol kg-1], ρ – the bulk
density of the soil [Mg m-3], D – the dispersion
coefficient [m2 s-1], and q – the water flux [m s-1].
The dispersion coefficient was calculated usingθ w D = λ q + θ w Dmτ w ,(2)

where λ is the dispersivity [m], Dm – the molecular
diffusion coefficient in free water [m2 d-1], and τw isa dimensionless tortuosity factor (τw = θw7/3/θs2)(Simunek et al., 1996b). The dispersivity and molecular diffusion coefficient are given in Tab. 1.

For the MIM simulations the CDE was modified
to consider two different types of soil water, a mobile fraction and an immobile fraction (vanGenuchten and Wierenga, 1976).

θ w = θ m + θim ,(3)

where θm is the mobile water content and θim – theimmobile water content. In our implementation of
the MIM model, the immobile water content during
the transient simulations was defined asθim = θ r + η (θ w − θ r ) ,(4)

where η is a dimensionless parameter varying between zero and one. We assumed that the residual
water (θr) was always immobile. The transport
equation for the MIM model is∂⎡θ m ck ,m + θim ck ,im ⎤⎦ +∂t ⎣∂
ρ ⎡⎣cˆk ,m + ck ,m + cˆk ,im + ck ,im ⎤⎦ =∂t
∂ck ,m⎤∂ ⎡− qm ck ,m ⎥ ,⎢θ w D
∂z ⎣
∂z
⎦

(5)

where the variables and constants have the same
meaning as in Eq. (1) but apply to either mobile(subscript m) or immobile (subscript im) water.
For the purpose of calculating solute transport,
the water flow velocities for the mobile water fraction were increased from values calculated with theRichards equation by a factor equal to θw/θm (qm =qwθw/θm). The MIM modification of the CDE involves the simultaneous solution of Eq. (5) for themobile fraction and a mass transfer equation expressing the transfer rate of solutes between the twowater fractions:T a b l e 1. Miscellaneous model parameter values.
T a b u ľ k a 1. Parametre modelu.
Parameter
1. Molecular diffusion coefficient in free water*
2. Longitudinal dispersivity**
3. Critical value for product of grid Peclet and Courant numbers***
4. CO2 concentration at soil surface
5. CO2 concentration at 3 m depth
6. Diurnal, sinusoidal temperature amplitude at soil surface
7. Maximum depth of water layer at soil surface

Value
1.0 x 10-4
0.05
2.0
3.6 x 10–4
0.02
5
0.05

Unit
[m2 d-1][m][–][PPMV][PPMV][°C][m]*Flury et al., 1999; **van Hoorn, 1981; ***Simunek et al., 1996b
1 – súčinitel molekulárnej difúzie v čistej vode, 2 – pozdĺžna disperzivita,3 – kritická hodnota súčinu Pécletovho a Courantovhočísla v uzle, 4 – koncentrácia CO2 na povrchu pôdy, 5 – koncentrácia CO2 v hĺbke 3 m, 6 – denná sinusoidálna amplitúda teploty na
povrchu pôdy, 7 – maximálna hrúbka vrstvy vody na povrchu pôdy.

77

P. J. Vaughan, P.J. Shouse, D.L. Suarez, M. Th. van Genuchten and J.E. Ayars∂θim ck ,im∂t= ω (ck ,m − ck ,im ) .

(6)

The MIM model has two additional parameters,
the immobile volumetric water content (θim) and themass transfer coefficient (ω) representing material
properties. Tab. 1 provides material property values
and other data that were kept constant in the simulations.
Eqs (5) and (6) were solved numerically bymeans of a one-dimensional Galerkin finite-element
approach similar to that described by Simunek et al.
(1998). The vertical flow domain at each location
was 3 m in length and was discretized into 377
elements. Nodal spacings varied exponentially from
3.0 x 10-3 m at the top to 0.03 m at a depth of 2.7 m.
Nodal spacings below this depth decreased linearlyin order to maintain stability in the transport calculations near the base of the simulated profile.
The mass balance error due to numerical approximations such as spatial and temporal discretizations can be calculated as the difference between
a true value, the sum of the integrated Cl- fluxes atthe boundaries of the column and the initial mass,and the current Cl- mass in the column. The maximum absolute value of the mass balance error for
the CDE model was found to be negligible. For all
MIM simulations the maximum error was 0.07
mmolc when the total mass in the column was 7.1
mmolc, for a relative error of 1%.
2.5 Field Crops and Irrigation
The field was cropped to tomato (Lycopersiconlycopersicum L.) in 1996 and cotton (Gossypium
hirsutum L.) in 1997. Root water uptake was simulated using
S (h, hϕ ) = α s (h)αϕ (hϕ ) S p ,(7)

where hϕ is the osmotic head [m], Sp – the potential
root water uptake distribution [m3 m-3 s-1] (Simunek
et al., 1996b), αs(h) – the water stress response
function, and αϕ(hϕ) – the osmotic stress response
function. The water stress response function was
calculated using

α s (h) = [1 + (h / h50 )b ]−1 ,(8)where h50 and b are crop-specific parameters(Simunek et al., 1996b). For tomato h50 = -5 m, b =2 and for cotton h50 =-25 m, b = 2. The values forcotton were estimated from results for alfalfa
78(Homaee et al., 2002). The value of h50 for tomato
was estimated by assuming a relatively high potential evapotranspiration rate of 6.4 mm d-1 (Hillel,
1998; p. 568). For osmotic stress, a similar equation, with the pressure head replaced by the osmotic
head (hφ), has the parameters hϕ50 = –100 m and b =3 for tomato and hϕ50 = –178 m and b = 3 for cotton(van Genuchten and Gupta, 1993). Rooting depth
was calculated as the product of the Verhulst-Pearllogistic growth function and the maximum rooting
depth (Simunek et al., 1996b). The maximum rooting depth for tomato was 1.2 m for tomato and 3.0
m for cotton. However, the program was modified
to reset the maximum rooting depth for cotton to
0.1 m below the water table.
Tomato was irrigated using gated-pipe water delivery to furrows, while cotton was furrow-irrigatedwith siphon-tube delivery at the upper end of the
furrow. Sprinklers were used for the pre-irrigations
in February of each year, prior to planting. The
precision of water delivery to the 65 ha field, as
determined from the precision of the gauge (0.1
acre-feet), was approximately 0.2 mm depth. Distribution non-uniformity was probably a greater
source of error in calculations of the infiltration
rate.
2.6 Initial Conditions
Initial volumetric soil water contents were calculated from measurements of the gravimetric watercontent and bulk density (Fig. 3). The gravimetricwater contents were determined by drying at 100°C
until no further weight change occurred. Initialpressure heads were calculated from water retentionfunctions utilizing hydraulic properties determinedwith the Rosetta program (Schaap et al., 1998).
Initial resident Cl- profiles were determined byanalysis of soil samples taken in November 1995
(Fig. 3).
2.7 Boundary Conditions
The atmospheric boundary conditions for waterflow at the soil surface alternatively involved timedependent prescribed fluxes and, if ponding occurred, prescribed pressure heads (Simunek et al.,1996b). Precipitation data were combined withirrigation data to build a continuous record of thepotential water flux at the soil surface for the twoyear simulation period starting November 1, 1995.
Calculations of the potential water flux assumed

Analysis of field-scale salinization: A comparison of different modeling approaches

Fig. 3. Measured initial resident Cl-(+) and volumetric water content (θw). Data shown are median values for 45 locations.
Obr. 3. Merané počiatočné rezidentné Cl-(+) a objemový obsah vody (θw). Údaje sú mediánom pre 45 lokalít.

that irrigation water was applied uniformly over thefield. While a necessary simplification for our onedimensional simulations, this boundary condition
likely was unrealistic for furrow-irrigation of both
crops. Possible development of ponding or runoff
was determined during the simulations. Irrigation inboth years was carefully managed to ensure thatwater did not run off the field. A maximum waterlayer of up to 0.05 m was allowed to build up,
thereby ensuring that no runoff was modeled. The
total depth of water applied to tomato, including the
pre-irrigation in February 1996, was 0.69 m. The
same total amount for cotton in 1997 was 0.85 m.
Transfer of water from the soil to the atmosphereoccurred by evaporation at the soil surface and byplant transpiration. Input data for calculating evaporation and transpiration rates were either daily potential soil evaporation and transpiration when thefield was cropped, or daily potential soil evaporation at other times. Both types of data were computed from CIMIS weather station data recorded at
the Panoche station some 20 km northwest of the
field. The daily potential evapotranspiration rate
was calculated as the product of a crop coefficientand the reference crop evapotranspiration rate reported each day by the station. Crop coefficients for
the San Joaquin Valley were taken from CIMIS
leaflet #214271). The CIMIS daily reference
evapotranspiration was calculated from a modified
Penman equation with a 0.10–0.15 m grass cover as
the reference crop (Doorenbos and Pruitt, 1977).
Potential soil evaporation rate was calculated fromthe reference crop evapotranspiration rate using an
expression for calculating bare soil potential evaporation from reference crop evapotranspiration (Snyder et al., 2000). The potential transpiration rate
was then estimated as the difference between thepotential evapotranspiration rate and the potential
soil evaporation rate.
The soil surface evaporation model requires
specification of a minimum value of pressure head(hcrit) at the soil surface (Simunek et al., 1996b). In
order to maintain numerical stability in the chemical speciation calculations we fixed hcrit at –25 m,
which at first caused too high estimates of the soil
1)

Available from DWR-DLA-WCO, 1020 9th Street, 3rd Floor,
Sacramento, CA 95814.

79

P. J. Vaughan, P.J. Shouse, D.L. Suarez, M. Th. van Genuchten and J.E. Ayars

evaporation rates. This problem was avoided by
specifying the potential soil evaporation rate after asoil wetting event using a square-root-of-time expression (Tanner and Jury, 1976; Black et al.,1969). During and for 2 days after an irrigation
event, potential soil evaporation rates were set
equal to the bare soil potential evaporation rate
(Snyder et al., 2000). After two days of drying, thepotential evaporation rate was reduced using thesquare-root-of-time expression with an empiricalconstant C = 0.4 (Black et al., 1969). Keeping hcrit
as –25 m during the first two days ensured that theactual soil evaporation rate matched the potential
soil evaporation rate at all times during the simulations.
The lower boundary of the simulated profile wasassumed to be 3 m below the soil surface, a depththat always remained below the fluctuating watertable. Pressure heads at this boundary were determined from groundwater levels in wells situated on
the same 6 x 6 grid as the sampling locations (Fig.
1). Weekly water table measurements were conducted during the growing season. At other times
the water table was estimated from fitted bicubic
surfaces using data from 18 permanent wells located within 3.3 km of the field. On three dates(July 1, 1996; July 9, 1997 and August 14, 1997)
we had simultaneous readings of the in-field wells
and the permanent wells. Mean water depths at welllocations in the field estimated from measurements
at the permanent wells were 1.50, 1.56 and 1.59 m
greater than the measured in-field values for thethree dates. Based on these results a correction of
1.55 m was applied to the water table estimatesfrom the bicubic surfaces. For water flow a
Dirichlet specified pressure head condition wasapplied to the lower boundary as determined fromthe water table measurements and estimates. For
chemical transport the lower boundary concentration was a specified concentration as determinedfrom chemical analyses of the groundwater at eachlocation. These concentrations showed little changeover the two-year period, thus allowing us to use aconstant concentration boundary condition.
3. Results and Discussion
3.1 Field Data Analysis

Resident chloride concentrations (molc m-3 soil)
for the final sampling in November, 1997 were
found to be highly variable among the 45 locations(for all measurements, CV = 142%). The arithmetic
mean of chloride concentrations of the six samples
80

taken from 0–1.8 m depth at each location varied
between 2.5 and 30.2 molc m-3. These mean values
were kriged using ordinary point kriging with aspherical, isotropic semivariogram model. The areal
distribution of the mean Cl- concentration exhibited
two maxima on either side of the northeastsouthwest trending band of relatively low concentration (Fig. 1). The lowest Cl- concentrations occurred along the eastern boundary of the field near
the main drain connecting the seven tile drains thatrun east-west. The low Cl- concentrations along theeastern boundary were likely influenced by somehorizontal transport to the sump in the northeast
corner.
3.1.1 Chloride Profile Shapes
Observed variations in the Cl- concentration withdepth may be classified by profile shape. Five different concentration distributions were identified:(1) monotonically increasing concentrations vs.
depth, (2) monotonically decreasing concentrations,distributions having (3) minima, or (4) maxima, and(5) other shapes, including sigmoidal shapes. The"maximum" shape is used to describe locations
where a maximum in the Cl- concentration occurredat some point within the top 1.5 m of the profile.
Figs 4 and 5 show examples of monotonically decreasing and maximum concentrations, respectively. The Cl- profile shapes were predominantly
monotonically increasing or of the maximum type(Tab. 2). Whereas the maximum shape was common for locations with high Cl- concentrations, thisshape was not found at any location where themean Cl- concentration for 0–1.8 m depth was lessthan 6 mmol L-1. The spatial distribution of the
maximum shapes was quite restricted. This shapeoccurred at all 9 locations in a triangular region
near the southeast corner delineated by a1, e5 andg1 in Fig. 1, and additionally in a region northwestof the center of the field, including locations g9, e9and j5.
The shapes of calculated Cl- profiles were similarly classified as above. Since the shape classification is subject to some uncertainty the classification
operation was repeated on three different days. The
largest variation in the shape classification data of
Tab. 2 was 10% for monotonically-increasing
shapes that were obtained with the MIM simulations.

Analysis of field-scale salinization: A comparison of different modeling approaches

3.1.2 Correlation of Cl- Concentration
Among Different Sampling Periods
Depth-averaged resident chloride concentrations
were calculated for each sampling except for thesecond sampling, in April, 1996, which was incomplete. Correlation coefficients were computed forthese data to characterize temporal variations. The
correlation coefficient matrix for the four complete

Fig. 4. Simulated in-situ Cl- concentration profiles [molc m-3] at location a3 for mobile (dotted) and immobile water (solid), andcomparison with measured concentrations at six depths recalculated in terms of the in-situ water content. Plot A is for the CDE, andplots B through I are for the MIM using MIM parameter values as shown.
Obr. 4. Rozdelenia koncentrácií Cl- [molc m-3] v poľných podmienkach v lokalite a3 pre mobilnú (bodkované) a imobilnú vodu
(plné), a porovnanie s meranými koncentráciami v šiestich hĺbkach, prepočítané na obsah vody in-situ. Tvar A je pre CDE a tvary B
až I sú pre MIM, použijúc MIM hodnoty parametrov, ako je to uvedené.

Fig. 5. Simulated in-situ Cl- concentration profiles [molc m-3] at location j5 for mobile (dotted) and immobile water (solid), andcomparison with measured concentrations. Plot A is for the CDE, and plots B through I are for the MIM using MIM parametervalues as shown.
Obr. 5. Rozdelenia koncentrácií Cl- [molc m-3] v poľných podmienkach v lokalite j5 pre mobilnú (bodkované) a imobilnú vodu
(plné), a porovnanie s meranými koncentráciami. Tvar A je pre CDE, a tvary B až I sú pre MIM, použijúc MIM hodnoty parametrov, ako je to uvedené.

81

P. J. Vaughan, P.J. Shouse, D.L. Suarez, M. Th. van Genuchten and J.E. Ayars
T a b l e 2. Cl- profile shape comparison [%].
T a b u ľ k a 2. Porovnanie tvarov profilov Cl-.
Monotonically increasing1)
Measured
49
CDE
91
MIM
44
1 – monotónne stúpajúce, 2 – monotónne klesajúce

Monotonically decreasing2)
400

samplings demonstrates that the time of year was
an important factor determining the correlation(Tab. 3). The highest correlation coefficient (r)were found for the first and fifth samplings, both inNovember. The other values of r that exceeded 0.7
occurred for post-harvest samplings. Lower valuesof r were generated when samples taken in thespring were compared with post-harvest samples.
This suggests that annual cyclic processes influenced Cl- concentrations in the vadose zone.
3.2 Simulations
3.2.1 CDE Modeling
Simulations were conducted for a two-year period starting on November 1, 1995.
Measured Cl- concentrations were determined for
0.3 m cores and, in order to compare the measured
value with model results, it was necessary to calculate mean modeled concentrations for depth intervals corresponding to those that were sampled.
Model results were compared with measured datafrom soil samples taken during four sampling periods starting in April 1996.
Initial simulations were conducted with potential
transpiration rates fixed as values estimated fromthe Panoche CIMIS station evapotranspiration data.
Considerable uncertainty exists in the CIMIS estimates of the reference crop evapotranspiration estimates since the CIMIS measurements were madeapproximately 20 km to the northwest. Calculated
resident Cl- concentrations for November, 1997,
using the CDE model with CIMIS potential transpiration rates, were smaller than measured concentrations at all depths and all locations. Furthermore, nomaxima were present in the Cl- distributions withdepth. This implied that upward transport of Clnear the soil surface was underpredicted. As a testof the effect of varying potential transpiration rates,
simulations were carried out using potential transpiration rates that were 1.5 and 2.0 times higher, but
without changing the potential soil surface evaporation rates.

82

Minimum
700Maximum
31
0
41

Other
9
9
15

Increasing the potential transpiration rate should
lead to less leaching because of more upward transport. However, an initial set of CDE results using
potential transpiration rates multiplied by a factorof 1.5 still produced overpredictions of Cl- leachingfor 38 out of 43 locations in November, 1997. The
five locations where observed leaching was greaterthan predicted were all near the eastern boundary of
the field (a7, a9, a11, ow-a, s). Simulations withdoubled potential transpiration did not exhibit significant improvement over those with a factor of
1.5. Hence, subsequent simulations discussed here
were performed with the potential transpiration
increased by 1.5.
We used the scaled root mean square error(SRMSE) to represent discrepancies between model
predictions and measured Cl- concentrations,
SRMSE =

100 ∑ ( Pi − Oi ) 2
,
N
O

(9)

where Pi are model predictions, Oi – observationsand ō is the arithmetic mean of N observations(Vanclooster et al., 2000). The median SRMSE for
the standard potential transpiration rate for 43 locations was 82%, while for otherwise identical conditions the mean SRMSE was 76% when the potential transpiration rate was increased by a factor of
1.5. The predictive capability of the CDE varied
substantially depending on location. One of thebetter predictions for all four sampling times wasfor location ow-s (Fig. 6). Nonetheless, even at thislocation, calculated resident Cl- concentrations in
November 1997 were significantly less than themeasured values at the four lowest sampling
depths.
The measured chloride concentrations exhibited
approximately lognormal distributions areallyacross the field at each of the six depths. A z-test
was conducted to determine the significance ofdifferences in the means of log10(Cl-) for the final
concentrations predicted with the model, and themeasured values from the November 1997 sampling. At α = 0.01 the test statistic exceeded thecritical value of 2.58, thus indicating a significant

Analysis of field-scale salinization: A comparison of different modeling approaches
T a b l e 3. Depth-averaged Cl- correlation coefficients.
T a b u ľ k a 3. Súčinitele korelácie Cl- spriemerované po hĺbke; 1 – dátum odberu vzoriek pôdy.
Sampling date1)
November, 1995
November, 1995
November, 1995
July, 1996
July, 1996
May, 1997

Sampling dateJuly, 1996
May, 1997
November, 1997
May, 1997
November, 1997
November, 1997r
0.73***
0.40**
0.87***
0.45**
0.78***
0.41**

Fig. 6. Profiles of the resident Cl- concentration [molc m-3] at location ow-s obtained with the CDE at four sampling times (as indicated in the lower left of the plots, where t is the number of days starting November 1, 1995). Circles are simulated concentrations
averaged over 0.3 m depth intervals. Measured data are shown by triangles.
Obr. 6. Rozdelenia koncentrácií Cl- [molc m-3] v poľných podmienkach v lokalite ow-s vypočítané s CDE pre štyri termíny odberov
vzoriek pôdy ( ako je to ukázané vľavo dolu na políčku, kde t je počet dní, počínajúc od 1. novembra 1995). Krúžkami sú označené
simulované koncentrácie spriemerované pre intervaly hĺbok 0,3 m. Merané hodnoty sú označené trojuholníkmi.

difference in the means for samplings in the 0–1.5
m range. For the lowest depth range, 1.5–1.8 m, thetest produced z = 2.19, which indicates that the null
hypothesis that the two means represented the same
population could not be rejected at α = 0.01. At thispoint the mobile-immobile model (MIM) features
were incorporated in the model to determine
whether the separation of Cl- between the mobile
and immobile water fractions would lead to better
simulations of the Cl- field data.
3.2.2 Mobile-immobile water modeling
Model parameters required for the MIM model
are the mass transfer coefficient (ω) and the immobile water content parameter (η). The tracer methods for obtaining MIM parameter values using tension infiltrometers (Clothier et al., 1992; Jaynes etal., 1995; Al-Jabri, 2002) could, at least in theory,
be applied to this field. However, given the large
spatial variations in texture, this would have required an inordinate number of measurements, wellbeyond the scope of our field measurement program. Instead, we elected to perform exploratory
simulations in which the MIM parameters werevaried systematically within certain ranges, i.e.,
10-4≤ ω ≤10-2 d-1, and 0.0 ≤ η ≤ 0.25. Optimal values of the two parameters for these ranges were
determined by locating the minimum in the SRMSE(Eq. (8)). We found that simulations could not be
reliably performed for η values above 0.25 without
inducing numerical instabilities due to high effective velocities of the mobile water. This problem
was more acute at locations having a highly variable saturated hydraulic conductivity of the different layers.
83

P. J. Vaughan, P.J. Shouse, D.L. Suarez, M. Th. van Genuchten and J.E. Ayars

Calculated in-situ Cl- concentrations were plottedfor various (η, ω) pairs and compared with measured in-situ Cl- concentrations sampled in November, 1997 at location a3 (Fig. 4). The solid line in
subplot A is the Cl- profile for a single CDE simulation. For the remaining MIM subplots, η had valuesof 0.0, 0.15 and 0.25 that increased horizontally
across Fig. 4, while ω varied within the range 10-4 ≤ω ≤ 10-2. Higher values of η yielded higher effective water flow velocities in both directions andcaused increasingly larger variations in Cl- for all
values of ω (Fig. 4). Concentration distributions
versus depth of Cl- in the mobile fraction were generally smooth, and diverged from Cl- in the immobile fraction as ω decreased. As to be expected, thecurves for the mobile and immobile Cl- concentration were almost identical when ω was relatively
large (ω = 10-2), thus reflecting rapid exchange
between the two liquid regions. Plots G,H and I
showed the smallest ω values (= 10-4) and exhibitedthe largest separation in Cl- profiles between themobile and immobile regions.
At location j5, neither model was able to reproduce the peak concentration observed in November
1997. In particular, this was the case for all MIM
simulations performed with varying values of η andω (Fig. 5). Application of the MIM to location ow-s
produced better agreement between the measured
resident Cl- concentration and the model results for(η = 0.25, ω = 10-4) as compared to other parametervalues and the CDE results (Fig. 7).
The MIM simulations produced a higher percentage of Cl- concentrations greater than measured
in the 0–0.3 m depth range (45% compared to 16%
for CDE). For the MIM the median SRMSE for thisdepth range was 22% compared to 41% for the
CDE, thus indicating that the MIM model performed substantially better than the CDE in theuppermost 0.3 m.
The SRMSE for each location was calculated forvarying values of the two MIM parameters. Contouring this statistic was a convenient way of studying the combined, interactive effects of the twoparameters on the calculations, including comparisons with the data (Fig. 8). Based on plots of thistype we obtained optimal values of the MIM parameters for each location. For 33 of the 36 locations on the 6 x 6 grid, the best agreement between
model results and data were obtained for the following combinations: 52% for (ω = 10-4; η = 0.25),
3% for (ω = 3.2 x 10-3; η = 0.25), 9% for (ω = 10-4;
η = 0.0), and 36% for (3.2 x 10-3 ≤ ω ≤ 10-2; η =84= 0.0). These results indicate that the MIM provided better agreement for low values of the transfer coefficient and/or high values of η at 64% of the
33 locations on the grid. For the remaining locations (36%), no significant improvement over the
CDE was apparent. The distribution of η was foundto be bimodal with 55% of locations having theminimum SRMSE when η = 0.25 and the remaining 45% when η = 0. No significant correlation was
found between the optimal MIM parameter values
and soil texture or bulk density. Three of the 36
locations were not considered in the MIM analysis
above since high immobile water contents, in association with rapid flow, caused numerical inaccuracies in the solute transport computations.
3.2.3 Cl- Profile Shape Classification
The mean weighted Cl- concentration of the
combined irrigation and precipitation fluxes during
the simulation period was 1.9 mmolc L-1, and themean lower boundary Cl- concentration for all locations, 17.7 mmolc L-1. This implies that solutions of
the transport equation, generally, should show anincrease in Cl- concentration with depth. Computed
Cl- profiles were mostly consistent with this as 91%
of the profiles for the CDE simulation were monotonically increasing (Tab. 2). However, only 49%of the measured Cl- profiles were of this type. Noneof the CDE profiles showed a maximum although
31% of the measured profiles had this shape. The
shapes of Cl- profiles calculated using the MIM
compared much better with the measured shapes(Tab. 2). In particular, 41% of the MIM simulations
generated maxima, albeit with lower peakconcentrations than those in the measured profiles.
The MIM predictions qualitatively were also muchcloser to the field data, although still with lower
maximum concentrations as exemplified by the
results in Fig. 5 for location j5.
3.2.4 Spatial Variations
Using the optimal MIM parameters for each location, the mean final Cl- concentration over the 1.8
m deep vadose zone was calculated. The results
were kriged using ordinary point kriging with aspherical, isotropic semivariogram model. Cl- distributions were similar to those for the measured
values (Fig. 9). Both maps had two maxima in approximately the same area of the field, and minima
along both the eastern boundary and southwest

Analysis of field-scale salinization: A comparison of different modeling approaches

Fig. 7. Profiles of the resident Cl- concentration [molc m-3], including Cl- in both water fractions, at location ow-s for November
1997. Plot A is for the CDE, and plots B through I for the MIM using MIM parameter values as shown.
Obr. 7. Rozdelenia koncentrácií Cl- [molc m-3] v lokalite ow-s pre obe frakcie vody v novembri 1997. Tvar A je pre CDE, a tvary B
až I sú pre MIM, použijúc MIM hodnoty parametrov, ako je to uvedené.

Fig. 8. Contouring of SRMSE at location ow-s based on simulations for 30 (η,ω) pairs represented by × symbols. The smallest
SRMSE occurred for the pair (0.25, 10-4).
Obr. 8. Obrysy SRMSE v lokalite ow-s zo simulácií pre 30 (η,ω) párov, reprezentovaných symbolmi ×. Najmenšie SRMSE bolizískané pre dvojicu (0.25, 10-4).

85

P. J. Vaughan, P.J. Shouse, D.L. Suarez, M. Th. van Genuchten and J.E. Ayars

Fig. 9. MIM model results: kriged surface of the mean Cl- concentration of the 0.0–1.8 m depth interval (November, 1997). Underlined labels signify Cl- maxima in the profile.
Obr. 9. Výsledky modelu MIM: krigovaný povrch priemerných koncentrácií Cl- v intervale hĺbok 0.0–1.8 m (november 1997).
Podčiarknuté značky označujú maximá Cl- v profile.

corner. The most noticeable difference was the lackof a distinct SW-NE trending minimum in thesimulated data. Instead, a saddle with no distinct
trend occurred near locations e5 and g5. The continuing, high Cl- concentrations at g9, m and surrounding locations could not be simulated using
either the CDE or MIM.
Considering that neither the CDE nor the MIM
model was capable of reproducing the complexityof the observed Cl- profiles, it is of interest to study
the spatial distributions of the residuals. For all four
sampling events subsequent to the initial sampling,
the CDE model performed slightly better (SRMSE= 46.6%) than the optimal MIM model (SRMSE =49.6%) for the six locations along the easternboundary of the field. The median SRMSE for CDE
at all other locations was 78.4% compared to 75.9%
for the MIM. The Cl- profiles with maxima tended
to be geographically associated with high mean Clconcentrations in the profile. Both the CDE andMIM significantly underestimated the field data atthese locations (Figs. 5 and 9). This suggests that
improvements may be needed in the representationof the flow and transport processes at the field site,or the boundary conditions involved, in future efforts addressing the development of Cl- concentration maxima.

86

3.3 Factors Potentially Affecting
the Predictive Capability of Both Models
Water tables in the field during the experimentwere managed as explained in the Methods section,i.e., groundwater levels were controlled by valvesinstalled at the eastern ends of tile drain laterals
while the field was cropped. Some uncertainty existed about the implications of this water table
management on the lower boundary conditions. We
believe that the overall effect was a reduction inleaching. Using a Dirichlet lower boundary condition for both the pressure head and the Cl- concentration may have been more appropriate for locations near the eastern boundary of the field. These
locations intersected a coarse-textured layer at
about 1 m depth and could have drained easily toward the sump in the northeastern corner of the
field, even with valves on the tile-drain laterals
stopping normal flow. This may explain the lower
SRMSE for these locations and the better performance here of the CDE model.
The accuracy of the mass balance of the modeled
profiles and the discrepancies between the model
predictions and the measured Cl- concentrations
imply that the simulated net downward fluxes ofchloride through the lower boundary were greaterthan the measured Cl- concentrations in NovemberAnalysis of field-scale salinization: A comparison of different modeling approaches

1997 would indicate. The shapes of the Cl- profiles
were consistent with this conclusion. We found that
35% of observed profiles were either monotonicallydecreasing or showed maxima. These shapes musthave been generated and maintained by greater
upward movement of Cl- than was predicted, particularly with the CDE model.
Mass balance discrepancies may have beencaused by at least three factors: (1) neither model
adequately represented the physical processes at thefield site, (2) either the lower or upper boundary
condition for water flow may have been poorly
represented for at least some portion of the 2-year
period, and (3) the root water uptake term in the
Richards equation was incorrectly calculated forsome depth range and time period during the simulations. Model inadequacy relates to the inability torepresent certain processes that are known to occur.
For example, a one-dimensional model cannot represent the effects of a sloping capillary barrier that
may add horizontal components to the flow velocity. While generally it may be easy to list the various processes that potentially may have been poorly
represented in a particular model, it is far moredifficult to assess the relative importance of those
processes to overall field-scale transport.
Our study suggests that the lack of specific data
concerning water table fluctuations and groundwater composition during the field experiment may
have contributed to a relatively poor definition of
the lower boundary condition. The upper boundary
condition was also subject to some uncertainty
since potential evapotranspiration rates were determined for a CIMIS site located 20 km away fromour field. Another possible source of uncertainty is
the representation of the potential soil surface
evaporation rate by a square root of time expression(Black et al., 1969). This simplification could be
especially significant in the off-season when baresoil potential evaporation (Snyder, 2000) accounted
for most or all of the upward water flux at and near
the soil surface. For example, tillage events during
the off-season could have caused substantiallyhigher evaporative fluxes for the same reasons thatthese events cause increased CO2 fluxes from soil
to atmosphere (Reicosky, 1997). Finally, irrigation
events were modeled by means of a uniform layerof water over the entire field. Some lateral nonuniform infiltration likely occurred due to local
variations in the hydraulic conductivity. Also, nonuniform infiltration may have occurred along the
furrows due to variations in opportunity times notaccounted for by the model.

The sink term in the Richards equation accounting for root water uptake at a specified depth was
determined using Eq. (7) involving the potential
transpiration rate, the root water uptake distribution
in the profile, and rooting depth. The potential transpiration rate had the same uncertainties with respect to applicability of the CIMIS data as the potential soil evaporation rate. While errors in the root
water uptake distribution and rooting depth may
have been important in calculations of the Cl- profile shapes, they probably did not contribute significantly to the overprediction of leaching.
Future studies at this site should include local
measurements of the latent heat flux to provide
better constraints on the calculated evapotranspiration rates. Also, measurements of water table elevations throughout the year could provide data toevaluate whether more detailed representations of
the lower boundary condition for water flow arenecessary. The same is true for the solute transport
boundary conditions near the water table.
4. Conclusions

Salinization of soils and salt loadings to groundwater are significant potential problems for agricultural productivity on the west side of the San Joaquin Valley. We attempted to model Cl- transport in
the upper 3 m of the soil profile at 45 locations in a
65-ha field in the San Joaquin Valley using both thestandard CDE and alternative MIM options in theUnsatchem program. Soil sampling was done atthese locations at five different times starting inNovember, 1995 and ending in November, 1997.
The predictive capabilities of the models were
evaluated by calculating the SRMSE for the fourlatter sampling events. The CDE model predicted
substantially greater leaching than actually occurred. We addressed this problem by upgrading
the model to account for MIM solute transport, andby increasing the potential transpiration rate by afactor of 1.5. Increasing the potential transpiration
rate caused more upward chloride transport, thereby
reducing the median SRMSE for the CDE from
82% to 76%. The MIM model simulations wereconducted by systematically varying the two MIM
model parameters (i.e., the immobile water fraction,
η, and the mass transfer coefficient, ω). A set ofsimulations utilizing optimal values of these parameters improved predicted Cl- profiles by reducing the median SRMSE to 72%. For 64% of thelocations, mass transfer coefficients within the87

P. J. Vaughan, P.J. Shouse, D.L. Suarez, M. Th. van Genuchten and J.E. Ayars

range 10-4 ≤ ω ≤ 10-3 together with relatively large
immobile water contents, provided the best predictions.
The CDE and MIM models were further evaluated by comparing the shapes of predicted andmeasured Cl- profiles. We found that the CDEmodel had a strong tendency to predict monotonically increasing resident Cl- profiles vs. depth for
November 1997. At the same time, measured chloride profiles exhibited a variety of shapes, including
31% having a maximum in Cl- within the top 1.5 m.
The MIM predicted maxima at 41% of the locations. Thus, predictions of profile shapes with the
MIM model were substantially better than such
shape predictions using the CDE. However, a significant deficiency of the MIM model was its inability to generate the large amplitude Cl- maxima
that were observed at some locations.
Future model studies may need to address the issue of how Cl- maxima are generated and maintained in irrigated soil profiles including the originof the areal distributions of the Cl- maxima. This
study also showed the importance of having accurate estimates of the evapotranspiration rate. For
example, continuous measurement of theevapotranspiration rate in a field experiment using
an eddy covariance system could provide a muchneeded constraint. Further studies of the lower
boundary condition of tile-drained fields with awater table transient are also warranted.
Ackowledgements. We acknowledge work done byRichard Schoeneman and Richard Soppe of the
USDA Water Management Research Laboratory in
Fresno. We are grateful to David Cone of the
Broadview Water District for providing records ofirrigation, precipitation and well depths. Salinity
Laboratory staff who assisted in this study were
JoAn Fargerlund and Jack Jobes.
REFERENCES
AL-JABRI S.A., HORTON R., JAYNES D.B., GAUR A.,
2002: Field determination of soil hydraulic and chemicaltransport properties. Soil Sci., 167, 353–367.
BELITZ K., PHILLIPS S.P., 1995: Alternative to agriculturaldrains in California’s San Joaquin Valley: Results of a regional-scale hydrogeologic approach. Water Resour. Res.,31, 1845–1862.
BENSON S.M., WHITE, A.F., HALFMAN S., FLEXSER, S.,
ALAVI M., 1991: Groundwater contamination at the
Kesterson Reservoir. California 1. Hydrogeologic setting
and conservative solute transport. Water Resour. Res., 27,
1071–1084.

88

BLACK T.A., GARDNER W.R., THURTELL G.W., 1969:
The prediction of evaporation, drainage, and soil water storage for a bare soil. Soil Sci. Soc. Am. Proc., 33, 655–660.
CLOTHIER B.E., KIRKHAM M.B., MCLEAN J.E., 1992:
In situ measurement of the effective transport volume forsolute moving through soil. Soil Sci. Soc. Am. J., 56, 733–
–736.
COTLOVE E., TRANTHAM H.V., BOWMAN R.L., 1958:
An instrument and method for automatic, rapid, accurate,and sensitive tritration of chloride in biologic samples. J.
Lab. Clin. Med., 51, 461–468.
DEVEREL S.J., GALLANTHINE S.K., 1989: Relation ofsalinity and selenium in shallow groundwater to hydrologic
and geochemical processes, western San Joaquin Valley,
California. J. Hydrol., 109, 125–149.
DEVEREL S.J. MILLARD S.P., 1988: Distribution and mobility of selenium and other trace elements in shallow groundwater of the western San Joaquin Valley, California. Environ. Sci. Technol., 22, 697–702.
DOORENBOS J., PRUITT W.O., 1977:Guidelines for predicting crop water requirements. Irrigation and Drainage Paper
24, 2nd Ed. Food and Agriculture Organization of theUnited Nations, Rome, 156 p.
DUST M., BARAN N., ERRERA G., HUTSON J.L., MOUVET C., SCHAFER H., VEREECKEN H. and WALKER
A., 2000: Simulation of water and solute transport in fieldsoils with the LEACHP model. Agric. Water Manage, 44,
1-3, 225–245.
FLURY M., YATES M.V. and JURY W.A., 1999: Numericalanalysis of the effect of the lower boundary condition on
solute transport in lysimeters. Soil Sci. Soc. Am. J., 63,
1493–1499.
HILLEL D., 1998: Environmental Soil Physics. Academic
Press, San Diego, CA.
HOMAEE M., FEDDES R A. and DIRKSEN C., 2002: Simulation of root water uptake II. Non-uniform transient water
stress using different reduction functions. Agric. Water
Manage., 57, 111–126.
JACQUES D., SIMUNEK J., TIMMERMAN A. and FEYEN
J., 2002: Calibration of Richards' and convection-dispersion
equations to field-scale water flow and solute transport under rainfall conditions. J. Hydrol., 259, 15–31.
JARVIS N., LARSSON M., 2001: Modeling macropore flow
in soils: Field validation and use for management purposes.
In: Conceptual Models of Flow and Transport in the Fractured Vadose Zone. pp. 189–215. National Academy Press.
Washington D.C.
JAYNES D.B., LOGSDON S.D., HORTON R., 1995: Field
method for measuring mobile/immobile water content andsolute transfer rate. Soil Sci. Soc. Am. J., 59, 352–356.
JURY W.A. and ROTH K., 1990: Transfer Functions andSolute Movement through Soil: Theory and Applications.
Birkhauser-Verlag. Basel, 226 pp.
MOHANTY B.P., BOWMAN R.S., HENDRICKX J.M.H.,
SIMUNEK J. and VAN GENUCHTEN M. T., 1998: Preferential transport of nitrate to a tile drain in an intermittentflood-irrigated field: Model development and experimental
evaluation. Water Resour. Res., 34, 1061–1076.
NIELSEN D.R., VAN GENUCHTEN M. T., and BIGGAR
J.W., 1986: Water flow and transport processes in the unsaturated zone. Water Resour. Res., 22, 89S–108S.
NKEDI-KIZZA P., BIGGAR J.W., SELIM H.M., VAN
GENUCHTEN M.T., WIERENGA P.J., DAVIDSON J.M.
and NIELSEN D.R., 1984: On the equivalence of two con-

Analysis of field-scale salinization: A comparison of different modeling approaches
ceptual models for describing ion exchange during transport
through an aggregated Oxisol. Soil Sci. Soc. Am. J., 46,
471–476.
REICOSKY D.C., 1997: Tillage-induced CO2 emission fromsoil. Nutrient Cycling in Agroecosystems. 49, 273–285.
SCHAAP M.G., LEIJ F.J. and VAN GENUCHTEN M.T.,
1998: Neural network analysis for hierarchical prediction ofsoil water retention and saturated hydraulic conductivity.
Soil Sci. Soc. Am. J., 62, 847–855.
SIMUNEK J., HUANG K., SEJNA M. and VAN GENUCHTEN M.T., 1998: The HYDRUS-1D software package forsimulating the one-dimensional movement of water, heat,and multiple solutes in variably-saturated media. Version
1.0, IGWMC-TPS-70, International Ground Water Modeling Center, Colorado School of Mines, Golden, Colorado,
186 pp.
SIMUNEK J., JARVIS N.J., VAN GENUCHTEN M.T. and
CARDENAS A., 2003: Review and comparison of models
for describing non-equilibrium and preferential flow andtransport in the vadose zone. J. Hydrol., 272, 14–35.
SIMUNEK J., SEJNA M. and VAN GENUCHTEN M.T.,
1996a: The HYDRUS-2D software package for simulatingwater flow and solute transport in two-dimensional variably
saturated media. Version 1.0, IGWMC-TPS-53, International Ground Water Modeling Center, Colorado School ofMines, Golden, CO, 167 p.
SIMUNEK J., SUAREZ D.L. and SEJNA M., 1996b: The
Unsatchem software package for simulating onedimensional variably-saturated water flow, heat transport,
carbon dioxide production and transport, and multicomponent solute transport with major ion equilibrium and kinetic
chemistry, Version 2.0. Research Report No. 141. US Salinity Laboratory, Riverside, California. 186 p.
SNOW V.O., CLOTHIER B.E., SCOTTER D.R. and WHITE
R.E., 1994: Solute transport in a layered field soil – experiments and modeling using the convection-dispersion approach. J. Contaminant Hydrol., 16, 4, 339–358.
SNYDER R.L., BALI K., VENTURA F. and GOMEZMACPHERSON H., 2000: Estimating evaporation frombare or nearly bare soil. J. Irrig. and Drain. Engng., 126,
399–403.
SUAREZ D.L., 2001: Sodic soil reclamation: Modelling andfield study. Aust. J. Soil Res., 39, 6, 1225–1246.
TANNER C.B. and JURY W.A., 1976: Estimating evaporationand transpiration from a row crop during incomplete cover.
Agron. J., 68, 239–243.
VAN GENUCHTEN M. Th., 1980: A closed-form equation forpredicting the hydraulic conductivity of unsaturated soils.
Soil Sci. Soc. Am. J., 44, 5, 892–898.
VAN GENUCHTEN M. Th. and GUPTA S.K., 1993: A reassessment of the crop tolerance response function. J. Indian
Soc. Soil Sci., 41, 730–737.
VAN GENUCHTEN M. Th. and SIMUNEK J., 2004: Integrated modeling of vadose-zone flow and transport processes. In: R.A. Feddes, G. H. de Rooij and J. C. van Dam(eds.), Unsaturated-Zone Modeling; Progress, Challenges
and Applications, Wageningen UR Frontis Series, Vol. 6,
Chapter 2, pp. 37–69, x-xi, Kluwer Academic Publishers,
Dordrecht, The Netherlands.
VAN GENUCHTEN M. Th. and WIERENGA P.J., 1976:
Mass transfer studies in sorbing porous media 1. Analyticalsolutions. Soil Sci. Soc. Am. J., 40, 473–480.
VAUGHAN P.J., SHOUSE P.J., GOLDBERG S., SUAREZ
D.J. and AYARS J.E.., 2004: Boron transport within an ag-

ricultural field: Uniform flow versus mobile-immobile water model simulations. Soil Sci., 169, 6, 401–412.
VAN HOORN J.W., 1981: Salt movement in unsaturated soil
profiles. In: A. Verruijt and F.B.J.Barends (eds.). Flow andTransport in Porous Media. Rotterdam, the Netherlands:
Balkema, 1981, 227–232 Proc. Euromech., 143. 2–4 Sept.
1981. Delft.
VANCLOOSTER M., BOESTEN J.J.T.I., TREVISAN M.,
BROWN C.D., CAPRI E., EKLO O.M., GOTTESBUREN
B., GOUY V. and VAN DER LINDEN A.M.A., 2000:
A European test of pesticide-leaching models: methodology
and major recommendations. Agric. Water Mgmt., 44,
1–19.
WAGENET R.J. and HUTSON J.L., 1986: Predicting the fate
and of nonvolatile pesticides in the unsaturated zone. J. Environ. Qual., 15, 315–322.
WAGENET R.J. and RAO B.K.,1983: Description of nitrogen
movement in the presence of spatially variable soil hydraulic properties. Agric. Water Mgmt., 6, 227–242.
WU L., SKAGGS T.H., SHOUSE P.J. and AYARS J.E., 2001:
State space analysis of soil water and salinity regimes in aloam soil underlain by shallow groundwater. Soil Sci. Soc.
Am. J., 65, 1065–1074.
Received 19. August 2004
Scientific paper accepted 8. February 2005

ANALÝZA ZASOLENIA V MIERKE POĽA:POROVNANIE ROZDIELNYCHMODELOVÝCH PRÍSTUPOV
P. J. Vaughan, P.J. Shouse, D.L. Suarez,
M. Th. van Genuchten, J.E. Ayars
Salinizácia pôdy a zaťaženie podzemných vôd soľami
je významným potenciálnym problémom poľnohospodárskej produkcie západnej časti San Joaquin Valley, California. Pokúsili sme sa modelovať transport Clvo vrchnej, 3m vrstve pôdy v 45 lokalitách na 65ha polív San Joaquin Valley, s využitím štandardných možností
CDE a alternatívnych MIM opcií programu Unsatchem.
Odber vzoriek pôdy v tejto lokalite bol vykonaný
v piatich termínoch, počnúc novembrom 1995 a končiac
novembrom 1997. Predpovedné možnosti modelu boli
vyhodnotené výpočtom SRMSE pre štyri po sebe nasledujúce termíny odberu vzoriek pôdy. Model CDE predpovedal podstatne vyššie vyplavovanie solí, ako bolo
v skutočnosti. Tento problém sme riešili zdokonalením
modelu pre MIM transport roztokov a zvýšili sme potenciánu intenzitu transpirácie 1,5-krát. Zvýšenie intenzity
potenciálnej transpirácie spôsobilo zvýšenie transportu
chloridov smerom nahor, čo spôsobilo zníženie mediánu
SRMSE pre CDE z 82 na 76 %. Simulácie modelom
MIM boli vykonané systematickou variáciou dvoch
parametrov modelu MIM (t.j. pomerným podielom imobilnej vody η a súčiniteľom prenosu hmoty ω). Séria
simulácií s využitím optimálnych hodnôt týchto parametrov zlepšila predpoveď profilov Cl- redukciou mediánu
SRMSE na 72 %. Pre 64 % lokalít koeficienty prenosu
hmoty v rozsahu 10-4 ≤ ω ≤ 10-3 spolu s relatívne vy-

89

P. J. Vaughan, P.J. Shouse, D.L. Suarez, M. Th. van Genuchten and J.E. Ayars

sokým podielom imobilnej vody viedli k našim najlepším prognózam.
Modely CDE a MIM boli potom vyhodnotené porovnávaním tvarov predpovedaných a meraných profilovCl-. Zistilo sa, že model CDE mal silnú tendenciu predpovedať monotónne klesajúci profil Cl- pre november
1997. V tom istom čase merané profily chloridov malirôzne tvary, 31 % z nich malo maximum Cl- vo vrchných
1,5m vrstvách. MIM predpovedal maximá v 41 % lokalít. Predpovede tvarov profilov s modelom MIM boli
podstatne lepšie ako predpovede s použitím CDE. Podstatným nedostatkom modelov MIM bola neschopnosť
generovať veľké amplitúdy maxím Cl-, ktoré boli pozorované v niektorých lokalitách.
Budúce štúdie modelov môžu byť nasmerované na
spôsob, ako sú generované a v závlahových podmienkach zachovávané maximá Cl-, vrátane dôvodu,
prečo tieto maximá Cl- sú priestorovo rozložené. Tátoštúdia tiež ukázala, aké je dôležité presné určenie intenzity evapotranspirácie. Napríklad, kontinuálne meranie
intenzity evapotranspirácie počas poľného experimentu
s využitím pulzačnej (eddy covariance) metódy môževiesť k veľmi potrebnému zlepšeniu výsledkov. Potrebný je tiež ďalší výskum dolných okrajových podmienok v pôdach so systematickou rúrkovou drenážou
v podmienkach neustáleného pohybu hladiny vody.
Zoznam symbolov
b
ck
ck
cˆk

90

– exponent vo funkcii vodného stresu [–],– koncentrácia k-tej zložky v pôdnom roztoku [M L-3],– koncentrácia k-tej zložky v substráte na povrchu[M M-1],– koncentrácia k-tej zložky v pevnej fáze [M M-1],

CV
D
Dm
h
hcrit

– koeficient variácie [–],– disperzný koeficient [L2 T-1],– koeficient molekulárnej difúzie vo vode [L2 T-1],– tlaková výška vody [L],– minimálna hodnota tlakovej výšky vody na povrchu
pôdy [L],hφ– osmotická tlaková výška [L],h50
– tlaková výška vody pri 50% redukcii stresom [L],hφ50
– osmotická tlaková výška pri 50% redukcii stresom[L],– pozorovaná koncentrácia Cl- [M L-3],
Oi– modelom predpovedaná koncentrácia Cl- [M L-3],
Piq– intenzita toku vody [L T-1],– potenciálna intenzita rozdelenia odberu vody koreňmiSp[L3 L-3 T-1],
SRMSE – relatívna hodnota kvadratickej odchýlky [–],t– čas [T],z
– vertikálna súradnica [L],αs
– funkcia závislá na vodnom strese [–],
αφ– funkcia závislá na osmotickom strese [–],
η– koeficient pre pre podiel “nepohyblivej” vody v pôde[–],θw
– objemová vlhkosť pôdy [–],
θr
– reziduálny obsah vody [–],
θs
– vlhkosť vodou nasýtenej pôdy [–],
θm
– vlhkosť „pohyblivej” časti vody v pôde [–],
θim– vlhkosť „nepohyblivej” časti vody v pôde [–],
λ– disperzivita [L],ρ− objemová hmotnosť pôdy [M L-3],
η– relatívna vlhkosť „nepohyblivej” vody v pôde [–],
τw
– tortuozita [–],
ω– koeficient prenosu hmoty [T-1].