Original Research

A. Yakirevich
T. J. Gish*
J. Šimůnek
M. Th. van Genuchten
Y. A. Pachepsky
T. J. Nicholson
R. E. Cady

Poten al Impact of a Seepage
Face on Solute Transport to a
Pumping Well
To develop predic ve models for solute transport in the subsurface, the mechanisms governing transport must be understood and quan fied. In this study, we used the HYDRUS-2D
variably saturated flow and transport model to describe the observed mass flux breakthrough curves (BTCs) of three surface-applied tracers—Br−, Cl−, and pentafluorobenzoic
acid (PFBA)—to a single well from which the tracers were pumped. Axisymmetrical transport simula ons of the data indicated the presence of an ac ve seepage face along thesoil–well interface near the water table. The calculated cumula ve water flux to the wellthrough the seepage face was found to be 120% of the variable flux boundary into the wellthrough the submerged zone. In addi on, calculated mass fluxes of Br, Cl, and PFBA throughthe seepage interface were approximately 8, 4, and 11 mes, respec vely, those throughthe variable flux boundary. Calcula ons suggest that a seepage face may be responsible forcausing the early arrival of the solutes and the overall shape of the BTCs. Our study indicates the poten ally important role a seepage face along the soil–well interface may havein modeling water flow and chemical transport to wells in variably saturated, unconfined
aquifers.
Abbrevia ons: BTC, breakthrough curve; MCP, moisture capacitance probe; PFBA, pentafluorobenzoic acid.

Simula on of groundwater flow and solute migra on are important forunderstanding the mechanisms governing chemical transport through the vadose zone
as well as lateral transport through the phreatic zone to surface water bodies. Modeling,
together with accurate monitoring protocols, may provide an integrated approach forassessing field-scale water and chemical fluxes through the unsaturated and saturated zones(Gish and Kung, 2007; Gerke et al., 2007).
Simula ng flux breakthrough curvesof three surface-applied tracers to asingle pumping well suggested thata seepage face along the soil–wellinterface in a variably saturated,unconfined aquifer can result in earlyarrival of the solutes.
A. Yakirevich, Zuckerberg Inst. for Water Science and Technology, J. Blaustein Inst. forDesert Research, Ben Gurion Univ. of the
Negev, Sede Boqer Campus, 84990, Israel; T.J.
Gish, USDA-ARS Hydrology and Remote Sensing Lab., Beltsville, MD 20705; J. Šimůnek,Dep. of Environmental Science, Univ. of
California, Riverside, CA 92521; M.Th. van
Genuchten, Dep. of Mechanical Engineering,
COPPE/LTTC, Federal Univ. of Rio de Janeiro,
UFRJ, Rio de Janeiro, 21945-970 RJ, Brazil;
Y.A. Pachepsky, USDA-ARS Environmental
Microbial and Food Safety Lab., Beltsville,
MD 20705; and T.J. Nicholson and R.E. Cady,
U.S. Nuclear Regulatory Commission, US NRCORR, Rockville, MD. *Corresponding author(Timothy.Gish@ARS.USDA.GOV).
Vadose Zone J. 9:686–696
doi:10.2136/vzj2009.0054
Received 5 May 2009.
Published online 3 Aug. 2010.

Different layouts may be used for tracer tests in the saturated zone (e.g., Domenico and
Schwartz, 1998, p. 506), including convergent radial flow and tracer transport toward a
pumping well. Unfortunately, such tests generally do not account for transport through theunsaturated zone and ignore the presence of a seepage face boundary, which is a common
feature of unconfined aquifers. A seepage face boundary is an interface at atmosphericpressure where the water exits the porous medium (Bear, 1972). In its vicinity, the phreatic surface is above the elevation of the external surface water level (e.g., an open drain, alake, or a pumping well). While seepage faces are often ignored in simulations of flow andtransport, they are believed to be important when the thickness of the unsaturated zone
is comparable to that of the unconfi ned aquifer (e.g., Shamsai and Narasimhan, 1991).
Several experimental and theoretical studies have highlighted the importance of seepage faces in controlling the vertical distribution of water and solute fluxes at or near the
groundwater–surface water boundary (Vachaud and Vauclin, 1975; Simpson et al., 2003;
Rushton, 2006; Chenaf and Chapuis, 2007). For example, in a sand box experiment ofradial flow and transport, Simpson et al. (2003) found that the flow velocity throughthe seepage face was about 30% higher than at the base of the aquifer, while the solutetransport patterns were strongly influenced by the flow characteristic of the system. They
concluded that numerical models ignoring the presence of a seepage face may incorrectlypredict solute residence time distribution.

All rights reserved. No part of this periodical may
be reproduced or transmi ed in any form or by anymeans, electronic or mechanical, including photocopying, recording, or any informa on storage andretrieval system, without permission in wri ng fromthe publisher.

To more accurately quantify solute fluxes through the unsaturated and saturated soil zones,
Kung et al. (2000) recently developed a method in which mass flux BTCs rather than soluteconcentration BTCs were collected from a tile drain. Th is tile drain protocol approachallows chemical transport processes to be precisely modeled during both unsaturated andsaturated flow conditions (Kladivko et al., 2001; Petersen et al., 2002; Gish et al., 2004;

2010, Vol. 9

www.VadoseZoneJournal.org | 686© Soil Science Society of America
5585 Guilford Rd., Madison, WI 53711 USA.

Köhne and Gerke, 2005). As a result, a large number of studies
have shown the potential for using mathematical models to simulate variably saturated flow and solute transport in tile-drainage
experiments (Villholth and Jensen, 1998; Mohanty et al., 1998;Šimůnek and de Vos, 1999; Larsson et al., 1999; Vogel et al., 2000,
Köhne et al., 2006; Gerke et al., 2007).
Gish and Kung (2007) recently extended the tile drain protocol toa non-tile-drained field where water and solute fluxes to a pumping
well through the unsaturated zone and a perched shallow aquifer
were monitored. Their procedure recovered >97% of three surfaceapplied tracers. Solute transit times based on a piston flow model
predicted that the tracers would reach the water table 4 d after
application, while observations showed the front arriving after
a little more than 2 d. Although preferential flow was initially
proposed as the reason for the early breakthrough times, anotherexplanation may be that the model did not account for a seepage
face at the soil–well interface near the water table. Since this well
protocol appears to be a useful tool for quantifying total solute
fluxes in non-tile-drained soils, it is important to compare resultswith more comprehensive numerical simulations that include
the effects of a seepage face. To our knowledge, most investigations assessing the impact of a seepage face on flow and transport
have been performed for steady-state flow conditions or without
accounting for water infi ltrating the vadose zone and recharging
the aquifer. The latter processes may well significantly affect prevailing flow patterns and solute fluxes toward a pumping well.
The objectives of this study were: (i) to analyze the Gish and Kung
(2007) tracer experiment using the HYDRUS-2D variably saturated flow and transport model (Šimůnek et al., 1999); (ii) toquantify water and solute fluxes through both the seepage zone
near the water table and the submerged portion of the pumping
well; and (iii) to evaluate a simplified analytical solution of the flow
equation for assessing the maximum application area for complete
tracer recovery during pumping (i.e., the well capture radius).

6 Materials and Methods

Tracer Experiment

Gish and Kung (2007) presented a detailed description of theirfield tracer experiment, which was performed at the USDA-ARS
Beltsville Agricultural Research Center in Beltsville, MD. Briefly,
the experiment was conducted in November of 2004 on a fieldthat had been under continuous corn (Zea mays L.) production
since 1996. The soil profile consisted of a coarse loamy sand surface
horizon (0–30 cm below ground level), followed by a sandy loam
horizon (30–80 cm), a coarse sandy horizon (80–150 cm), a gravelly sand layer (150–250 cm), and finally a sandy clay loam to clay
loam horizon (250–350 cm) of very low permeability. The average
depth of the perched water table at the site during the week before
the experiment was 168 cm.

Fig. 1. Sketch of the field experimental design (not to scale).

A schematic of the experimental design is presented in Fig. 1. The
plot was irrigated for 7 d at a constant rate of 4.1 mm h−1 using
seven rows of four solid-state sprinklers spaced 6.4 m apart. Water
pumping was applied using a gasoline pump installed in the observation well (5-cm diameter). The well was constructed with slotted
polyvinyl chloride from 0.5 to 3.0 m. Once the slotted well wasinserted and well extensions attached, blow sand was poured into
the gap between the well and the native soil (for just the 5-cm hole).
As slots in the well can become clogged during soil insertion, positive water pressure was applied to the top of the observation well
to clean out the well slots and activate the well. The water level was
continuously monitored with a submersible pressure transducer(Eijkelkamp Agrisearch Equipment, Giesbeek, the Netherlands)
located at 2.9 m below the soil surface in the observation (pumping) well. A line for pumping water from the observation well wasinserted 0.15 cm below the water table depth and attached to a
portable gas pump, which ran throughout the experiment. The
setup created a converging flow field in the saturated zone toward
the observation well. Moisture capacitance probes (MCPs) were
installed 11 m away from the pumping well to monitor soil watercontents every 10 min at depths of 0.1, 0.3, 0.5, 0.8, 1.2, 1.5, and
1.8 m below the soil surface. Three conservative tracer pulses were
uniformly manually sprayed between radii of 0.1 and 3.0 m aroundthe well, resulting in a treated area of about 28.2 m2 (Fig. 1): Br−,
Cl−, and PFBA. The Br pulse (0.1064 kg Br m−2) was applied asa surface broadcast spray 1.5 d before irrigation was initiated, Cl(0.052 kg Cl m−2) was applied as a surface broadcast spray immediately before the irrigation, and PFBA (0.0167 kg PFBA m−2)
was applied as a surface broadcast spray 7 h after the experiment
began. Each tracer was applied in just enough water to get it into
solution. The background pumping rate before irrigation was about
0.03 L s−1 and reached a maximum rate of 0.18 L s−1 at the endof the study. Immediately following Cl application, the irrigation
system was turned on and water samples were manually collected
every 15 to 60 min for the duration of the study. We used theobserved water content data to estimate the soil hydraulic parameters needed for the unsaturated flow model, while the tracer BTCs
in the pumping well were used to evaluate the performance of thesolute transport model.

www.VadoseZoneJournal.org | 687

Mathema cal Model
The tracer experimental data were analyzed using the
HYDRUS-2D model simulating variably saturated water flow andsolute transport (Šimůnek et al., 1999). The model is based on the
Richards equation for flow and the advection–dispersion equation
for transport. Assuming a conservative tracer (no sorption anddecay), these equations are, respectively:∂θ= ∇ ⎡⎣ K ( ∇ψ + 1)⎤⎦
∂t

[1]

∂θ C
= ∇ ( θ D ∇ C − qC )∂t[2]

where θ is the volumetric water content [L3 L−3], t is time [T], ψ
is the pressure head [L], K is the hydraulic conductivity tensor [LT−1], C is the solution concentration [M L−3], D is the hydrodynamic dispersion tensor [L2 T−1], and q is the Darcy–Buckingham
flux vector [L T−1].
The soil water retention and the unsaturated hydraulic conductivity functions needed for the solution of Eq. [1] were described using
the van Genuchten (1980) relationships:Se =θ − θrn −m= ⎢⎡1 + (α ψ ) ⎥⎤⎦θs − θr ⎣

[3]

m ⎤2
⎡
K ( Se ) = K s Seτ ⎢1 − (1 − Se1/m ) ⎥
⎣⎢
⎦⎥

[4]

where S e is relative saturation, θ r is the residual water content[L3 L−3], θ s is the water content at saturation [L3 L−3], Ks is thesaturated hydraulic conductivity [L T−1], α [L−1] and n are shapeparameters, τ = 0.5, and m = 1 − 1/n. The density effect was notaccounted for because measured concentrations of solutes did notexceed 1 g L−1.
Application of Eq. [1–4] to the field tracer experimentusing HYDRUS-2D requires definition of the numerical finite element mesh, the initial and boundary
conditions, and the soil hydraulic and solute transport parameters. We considered an axisymmetricalflow field of 250-cm height and a radius of 3000 cm(Fig. 2). Vertical extension of the simulation domain
was chosen due to the existence of a nearly impermeable clay layer at the depth of 250 cm that prevented
deep leaching and caused lateral water movement andsolute transport. The grid origin (radius r = 0, depth z= 0) was set at the center of the observation well andat a depth of 250 cm from ground level. The finite element mesh size was 5 cm in the vertical direction and10 cm in the horizontal direction everywhere, except near the soil
surface and the well boundary, where we used a mesh size of 2.5
cm. The low-permeability clayey bed at the bottom of the simulation domain was considered to be a no-flux boundary for bothflow and solute transport. Along the right boundary, 30 m fromthe pumping well, we assumed that the water table depth (1.68 m)
was not affected by the irrigation by imposing an equilibrium pressure head distribution vs. depth for flow and a zero concentration
gradient for transport.
The soil surface boundary was divided into four concentric rings(Fig. 2). A constant water flux of 4.1 mm h−1 was imposed between
radii of 2.5 to 2000 cm and a zero flux elsewhere. A solute flux withconcentration equal to the concentration of the irrigation water(Cir) was assigned between radii of 2.5 and 10 cm and between 300
and 2000 cm. Water and solute fluxes were assumed to be zero forradii >2000 cm. Flux concentrations between 10 and 300 cm near
the well where the tracers were applied were calculated such thatthey were consistent with the total tracer application rates. We
assumed pulse durations of 1 h (related to the dissolution times)
for all tracers, leading to concentrations of the infi ltrating waterof 25.95, 12.69, and 4.07 g L−1 for Br, Cl, and PFBA, respectively.
The Cl concentration of the irrigation water was 8.75 mg L−1 (Gishand Kung, 2007), while no Br or PFBA were present in the irrigation water.
Along the left boundary (at r = 2.5 cm) representing the well, a
seepage face was allowed to develop above the water level in thewell. Although the model stipulated a potential seepage boundaryfrom the elevation of 85 cm to the soil surface, the “active” seepage face occurs only near the water table (see Fig. 2), dependingon the transient flow conditions during the experiment. During
simulation, the saturated part (ψ ≥ 0) of a potential seepage faceis treated as a prescribed pressure head boundary with ψ = 0, while
the unsaturated part (ψ < 0) is treated as a prescribed zero flux
(Fig. 2). The lengths of the two surface segments are continually
adjusted (Neuman, 1975) during the iterative process until theFig. 2. Simulation domain with the invoked boundary conditions (not to scale). The
profile consisted of five different layers (Layer 2 is divided into two sublayers; Layer
5 is composed of low-permeability clay loam and is not included in the simulation
domain); CT and Cir are solute concentrations of the infiltrating water.

www.VadoseZoneJournal.org | 688

calculated values of flux along the saturated part, and the calculated values of ψ along the unsaturated part, are all negative, thus
indicating that water is leaving the flow region through the saturated part of the surface boundary only.
Additionally, a variable flux (Neumann boundary condition) was
specified between the elevations of 5 and 80 cm in the submerged
well zone (Fig. 2) where the pump was installed. A more rigorous
procedure would be to specify equilibrium pressure head distribution (Gardner, 2005; Wilson and Gardner, 2006) based ontransient water levels measured by the pressure transducer; however, the obtained data exhibited a highly oscillating water level
due to non-uniform pumping rates and relatively small water volumes in the well, which makes it difficult to prescribe transientpressure head values. Another reason why we used the variable flux
boundary instead of a prescribed pressure head is that we wanted topreserve water fluxes pumped from the well. Note that we did notneed to account for any loss in flows due to plugged well slots. The
imposed variable flux through the submerged well boundary was
adjusted iteratively during the numerical calculations. We started
by prescribing the flux values equal to the pumping rate, and calculated the total water volume out of the well boundary (a sum of theseepage flow and the variable flux), which significantly exceeded
the total volume pumped. We gradually decreased the imposedvalues of variable flux until the simulated total flow volume along
the well boundary matched the measured total volume of pumped
water at the end of the simulation period. Finally, a zero concentration gradient was assumed along the left and right boundaries of
the transport domain.
As the initial condition for the pressure head we used an equilibrium distribution in the profi le corresponding to the observed
initial water table elevation of 82 cm. The initial concentrations
for the transport simulations were zero for Br and PFBA and 7.3
mg L−1 for Cl, being the average measured Cl concentration in thesaturated zone before irrigation.

Calculated Tracer Concentra ons in the
Pumping Well
Solute concentrations of the pumped water (Cw) were calculated
from a mass balance equation assuming full mixing in the pumping well as follows:
π rw2 hw

dC w
= Q c(SF) (t ) + Q c(VB) (t ) −Q p (t )C w
dt

[5]

where r w is the well radius [L], hw is the time-averaged mean waterlevel in the well measured from the well base [L], Qc(SF) and Qc(VB)
are the mass fluxes [M T−1] of the tracers to the well through theseepage face and variable flux boundaries, respectively, as calculated with HYDRUS-2D, and Q p is the pumping rate [L3 T−1].
We used the time-averaged mean water level in the well because
the actual record exhibited highly oscillating values. In fact,

simulations indicated that the effect of hw on concentration was
not significant due to the relatively small water volume in the welland the large water and solute fluxes.
The goodness-of-fit of the model was evaluated using the coefficientof determination (R2) and the modified index of agreement (MIA)
as given by Legates and McCabe (1999):
MIA =n∑ i =1 C w (ti ) −C wobs (ti )1−nn∑ i =1 C w (ti ) −C wobs + ∑ i =1 C wobs (ti ) −C wobs

[6]

where Cw(ti) and Cwobs(ti) represent the simulated and observed
concentrations of a tracer in pumped water at time ti, and C wobs isthe mean observed concentration. While there is no statistical basis
to decide exactly which MIA value is a good threshold characterizing the use of an “accurate” model, following Köhne et al. (2005)
we assumed simulations with MIA > 0.75 as being “accurate.”

Inverse Solu on for the Unsaturated Soil
Hydraulic Parameters
The soil hydraulic parameters needed for the calculations were estimated using a combination of one- and two-dimensional inverse
simulations performed with the HYDRUS-1D (Šimůnek et al.,2005) and HYDRUS-2D soft ware packages, respectively, as well aslaboratory measurements of the saturated water content (θs) of thevarious soil horizons. During a first step, we used HYDRUS-1D toestimate the hydraulic parameters from the observed water content
values that were monitored at different depths in the soil profi le
11 m away from the pumping well. By using HYDRUS-1D weassumed that flow in the unsaturated zone relatively far (11 m)
from the pumping observation was mainly in the vertical direction, thus considerably minimizing computational times comparedwith using HYDRUS-2D. During the second step, we used
HYDRUS-2D to perform additional calibrations of the saturated
hydraulic conductivity (Ks) of Layers 3 to 5 (Fig. 2), where waterpresumably flowed both vertically and horizontally.
For the HYDRUS-1D simulations, we assumed a 250-cm-deepsoil profi le consisting of five layers consistent with the described
lithology of the site: loamy sand at 0–30 cm, sandy loam at 30–50
cm, sandy loam at 50–80 cm, coarse sand at 80–140 m, and gravely
sand at 140–250 m. The clay loam horizon (250–350 cm) had an
extremely low conductivity and as such was considered impermeable and hence not included in the simulations. The sandy loam
horizon (30–80 cm) was subdivided into two layers based on thefact that markedly different water contents were observed during
steady infi ltration at depths of 30 and 50 cm, thus indicating different hydraulic properties. The MCP-measured water content at
depths of 10, 30, 50, 80, 120, and 150 cm were used to calibrate theflow model, which in total contained 25 hydraulic parameters (thefive parameters θ r, θ s, α, n, and Ks for each of the five layers). To

www.VadoseZoneJournal.org | 689

minimize issues of uniqueness in the inverse solution, we decreased
the number of optimized parameters by using experimentallymeasured values of the water content at saturation (θ s), and zero
values of the residual water content (θ r) for all coarse-textured
layers (except sandy loam). Thus, we needed to determine threeparameters (Ks, α, and n) for each of the five layers and one θ r value
for the sandy loam horizon.
The experiment was conducted in November of 2004 when air
temperatures ranged from 1 to 7°C, the evaporation rate waslow, and no rain or irrigation occurred in the week before theexperiment (Gish and Kung, 2007). Therefore, we assumed thatthe initial pressure head profi le was close to equilibrium with thewater table at 168 cm. The boundary condition at the soil surface
defined a water flux of 4.1 mm h−1 during the irrigation period.
Since the bottom boundary at 250 cm was considered to be impermeable, the water table was likely to rise at some distance aroundthe well after irrigation began. As a result, some lateral water flowin the saturated zone occurred, which had to be accounted forby the flux through the lower boundary in the one-dimensional
HYDRUS-1D simulations. This flux varied with time, dependingon the position of the water table. To simulate this process, we usedthe HYDRUS-1D horizontal drain boundary condition with adrain spacing of 22 m (twice the distance from the MCP location
to the pumping well). Although this boundary condition does
not describe the physical two-dimensional flow system exactly, it
was the boundary condition available in HYDRUS-1D that mostclosely represented our field conditions.
Initial estimates of the α and n parameters for the inverse procedure were based on the assumed initial equilibrium distributionof the pressure head and the observed water contents before theexperiment, while Ks values were initially estimated from soil
textural class and particle size distribution information. We useda sequential inverse procedure by starting the parameter searchwith the first layer, while assuming that the other parameters were
known (the initial estimates) and using observed water contentsof the 10-cm depth only. The next step was to find parameters of
the second layer by fitting water contents at depths of 0.1 and 0.3m, and using the parameters of the first layer found at the previous step. The third step was to adjust the parameters of Layers 1
and 2 simultaneously by fitting the measured water contents at
depths of 0.1 and 0.3 m. We continued this sequential procedureby gradually increasing the number of layers and hence the numberof parameters that were simultaneously estimated. During the final
step, the code searched simultaneously for 15 parameters (θr for the
sandy loam horizon was not changed after the third run).

Well Capture Zone under
Steady-State Condi ons
The recovery of the tracers is an important factor that allows assessment of the solute mass balance. In a radially convergent, saturated
flow field, the recovery depends on the tracer application radius

and the well capture zone. To avoid losses of the tracer with lateral
flow, the radius must be less than the radius of the water tabledivide, where the radial flow velocity is equal to zero (Fig. 3). The
radius of the water table divide can be roughly estimated from asolution of the Boussinesq equation for steady-state radial unconfined flow toward a well, subject to Dirichlet boundary conditions
as follows:
1 ∂ ⎛⎜∂ h ⎞⎟⎟ = −iH ( a − r ) r0 ≤ a ≤ R⎜⎜⎝ K s rhr ∂r
∂ r ⎠⎟[7]r = r0 : h = h0

[8]

r = R : h = hR

[9]where h is the water table height [L], Ks is the hydraulic conductivity of the aquifer, i is the recharge rate [L T−1], H(r) is the Heaviside
step function, a is the radius of the infi ltration area (r 0 ≤ a ≤R), r 0 and R are the radii of the well and the outer boundary [L],respectively, h0 and hR are the water table heights at r = r 0 and r =R, respectively, and r is the radial coordinate.
Assuming a constant hydraulic conductivity, Ks, and integrating
Eq. [7], we obtained an analytical expression for the radial water
flux to the well, Q:
∂h
∂r
2
⎧
⎫
⎪
⎪
ir
H
a − r ) − ia 2 H ( r − a )
(
⎪
⎪
⎪
⎪
⎪
⎪
2
2
⎛
⎞
⎪
⎪
a
−
r
R
⎟
2
2
2
⎜
0
⎪
hR − h0 ) K + i ⎜⎜+ a ln ⎟⎟ ⎪
= π⎨
(
⎬
⎟
⎪
⎜⎝ 2
a ⎠⎪
⎪
⎪
⎪
⎪
−
⎪
⎪
⎪
⎪
ln ( R r 0 )
⎪
⎪
⎪
⎪
⎩
⎭

Q = − 2π Krh[10]Fig. 3. Schematic of unconfined radial flow with infiltration; r, radial
coordinate; a, radius of infiltration area; r*, well capture radius; h,water table height; r0 and R, radii of well and outer boundary, respectively; h0 and hR, water table heights at r = r0 and r = R, respectively;
Q, radial water flux.

www.VadoseZoneJournal.org | 690

For the limiting case of a = R, this solution coincides withan equation derived by Aravin and Numerov (1965).
An equation for the water table divide radius (r*) may beobtained by solving the equation for Q = 0 for the casewhen r* < a:
⎛ 2

r* =

⎞

( hR2 − h02 ) Ki s + ⎜⎜⎜⎜⎝ a −2 r0 + a 2 ln Ra ⎠⎟⎟⎟⎟
2ln ( R r0 )[11]

The pumping well discharge is calculated from
Q r0 = − 2π r0 K s h0

∂h
r =r0
∂r

⎛ a2 −r 2⎪⎧⎪R ⎞ ⎪⎫
0⎪⎪hR2 − h02 ) K s + i ⎜⎜⎜+ a 2 ln ⎟⎟⎟ ⎪⎪⎪(⎜a ⎠⎟ ⎪⎪⎪⎝ 2= π ⎪⎨ir02 −⎬ [12]⎪⎪⎪⎪ln ( R r0 )
⎪⎪
⎪⎪
⎪⎪
⎪⎪
⎩
⎭

Fig. 4. Observed (circles) and simulated (lines) water contents at different depths.

6 Results and Discussion

Unsaturated Soil Hydraulic Parameters

Figure 4 compares the MCP-measured water contents at different
depths with the simulated values obtained with HYDRUS-1D.
Good agreement was obtained between the observed and simulated water contents for all depths except at 0.8 m, for which R 2= 0.877. The fluctuations in observed water contents at the four
upper observation points may have been caused by non-uniform
irrigation during several windy days (Gish and Kung, 2007). The
calculated water content time series clearly describes the propagation of the moisture front into the profi le.
Table 1 presents the fitted hydraulic parameters obtained using the
invoked inverse procedure. The values of α and n are within therange reported by Hodnett and Tomasella (2002) for soils withsimilar textures. A relatively high value of 10 was obtained for thevan Genuchten hydraulic parameter n for the fi fth layer (gravelysand), which implies a very steep water retention curve and an
unsaturated hydraulic conductivity curve in which the conductivity decreases very rapidly with decreasing pressure head. Whilesome slight changes in the minimized objective function occurred
when the parameter n varied between 8 and 12, the simulated
water contents were found not to be sensitive to these changes.
We therefore accepted a value of n = 10 for the very coarse-textured
gravelly sand layer. The saturated hydraulic conductivity for the
loamy sand and sandy loam (Layers 1 and 2, Table 1) seemingly
has extremely high values, yet these values were obtained as a resultof inverse simulations for the limited range of water content in
the unsaturated zone and for the most common value of the poreconnectivity parameter τ = 0.5 (Eq. [4]). The value of τ can vary
in a range from −1 to 2.5 (Mualem, 1976), however, and Schuhand Cline (1990) reported that τ varied between −8.73 and 14.80.
This most significantly affects the hydraulic conductivity values atwater content close to saturation. Thus, the obtained Ks values for
Layers 1 and 2 have a great deal of uncertainty and should not be
used to simulate saturated conditions.
Additional calibrations were performed with the complete twodimensional axisymmetrical flow model against the observed watercontents. As a result, we obtained values for Ks of 185, 149, and 52
cm h−1 for Layers 2-2 (sandy loam), 3 (coarse sand), and 4 (gravelysand), respectively. The maximum difference of these Ks values
from those obtained for the one-dimensional problem was for the
fourth layer. The calibrated Ks for the two-dimensional problem
here was approximately one-third than that of the one-dimensional
problem (Table 1). This difference was probably due to using animperfect lower boundary condition in the one-dimensional case.
The simulated water contents obtained with HYDRUS-2D were
essentially identical to those shown in Fig. 4. For the two-dimensional calibration R 2 = 0.858, slightly smaller than the value of
0.877 obtained with the HYDRUS-1D calibration.

Table 1. Unsaturated soil hydraulic parameters† inversely estimatedusing HYDRUS-1D.
Layer

Depth

θ s‡

θr

α

n

cm−1

cm
0.327

0

0.046

Ks
cm h−1

1

0–30

1.62

100

2-1

30–50

0.351

0.024

0.054

1.46

115

2-2

50–80

0.260

0

0.054

1.33

113

3

80–150

0.315

0

0.069

1.76

161

4

150–250

0.330

0

0.044

10.00

167

† θ s and θ r, saturated and residual water contents, respectively; α and n, shapeparameters; Ks, saturated hydraulic conductivity.
‡ Experimentally determined.

www.VadoseZoneJournal.org | 691

Despite the achieved good fit to the experimental data, we realizethat the obtained set of parameters may not be unique and that
additional experiments would be helpful to reduce uncertainty
in the parameters. Nevertheless, we accepted these values andthen used them for our simulations for the two-dimensional case.
Of course, the objective of our study was to obtain an accurate
description of the experimental field data (water contents, flow
rates, and concentrations of the well water), rather than uniquesoil hydraulic parameters. This is reasonable since no attempts weremade to extrapolate the data observed beyond the experimental
time period.

Two-Dimensional Simula ons of Water Flowand Solute Transport
Solute transport simulations were performed for four different
values of the longitudinal dispersivity, aL (5, 15, 25, and 40 cm).
These values are within the range of 0.05 to 1 m obtained from thefield tracer tests (Gelhar et al., 1992). We used the same value forall layers, while assuming in most cases a ratio of 10 between thelongitudinal (aL) and transversal (aT) dispersivities. This ratio iswell within the range of ratios (5–20) reported in several studies byAnderson (1979) and Domenico and Schwartz (1998, p. 506). For
unsaturated soils, the variability in flow direction and microscopic
velocity can be larger than in saturated soils because of the lower
degree of water saturation (Maraqa et al., 1997; Nützmann et al.,2002). A clear relationship between the degree of saturation and thedispersivity has not yet been established, however, and therefore weused constant values of these parameters for the entire variably saturated domain. Figure 5 compares the observed and calculated (using
Eq. [5]) solute BTCs of the tracers in the well. The results are shownfor simulations with and without accounting for the seepage face.
Neglecting the seepage face boundary caused a delay of a few hours
in the arrival of the solute fronts as well as of the maximum concentration, except those of PFBA. The PFBA was applied 7 h after
the experiment began and the flow conditions in the unsaturated
zone were different than at the beginning. The peak concentrations
decreased somewhat, except perhaps for Cl (Fig. 5b). Because overall
agreement with the experimental data was better when the seepage
face was considered, we will focus primarily on results obtained for
simulations with the seepage face boundary.
The best agreement between the observed and simulated BTCsof the tracers in the well was obtained when we used a value of 15
cm for aL (Fig. 5). The R2 values were 0.909, 0.899, and 0.930 and
the MIA criterion (Eq. [6]) was 0.883, 0.857 and 0.896 for Br,
Cl, and PFBA, respectively. The model provided a good description of the tracer arrival in the well, but overestimated the arrival
time of the maximum concentration by around 8.7%. Simulations
with a dispersivity of 5 cm improved the description of the front
part of the solute BTCs (especially the initial slope) but caused
much lower tails. On the other hand, simulations with a dispersivity of 40 cm described the tails quite well but resulted in muchearlier arrival times and lower maximum concentrations than theFig. 5. Observed (circles) and simulated solute breakthrough curves of(a) Br, (b) Cl, and (c) pentafluorobenzoic acid (PFBA). Simulations
were obtained with (solid lines) and without (dashed lines) accounting for a seepage face boundary.

observed values. Decreasing the ratio aL/aT from 10 to 5 did nothave an appreciable impact on the simulated BTCs.
The steeper arrival fronts, more pronounced asymmetrical shape,earlier concentration peaks, and longer tails of the observed BTCs
compared with the simulated curves (Fig. 5) may indicate enhanced
preferential flow or an underestimation of the soil heterogeneity
in the simulations (Kung et al., 2000; Buczko and Gerke, 2006).

www.VadoseZoneJournal.org | 692

Some of these features probably could be accounted for by invoking a scale-dependent dispersivity in the model; however, limited
soil data were available. The variability in the observed concentrations at about the time when the peak concentrations were detected(Fig. 5) may have been caused by an unsteady pumping rate (thepumping rate generally increased with time) and because the pump
had to be replaced at some point during the experiment (Gish andKung, 2007). The simulated BTCs (Fig. 5) also exhibited oscillating behavior during the time interval when the pumping rate
fluctuated significantly. The latter caused changes in the seepage
face length and in the values of solute fluxes through the seepage
face and the submerged part of the well. Abrupt changes in thesefluxes and in their ratio induced the oscillation of simulated concentrations in the well.
The sensitivity of the model to the assumed dissolution times of the
tracers at the soil surface was also evaluated. A fourfold decrease or
increase in the pulse duration (0.25 or 4 h, respectively) and consequently a four times increase or decrease in the surface boundary
concentration (to preserve the injected tracer masses) were foundto produce similar results. We additionally tested the effect of
having a different thickness of the fi fth layer (where most of thewell screen was located) by moving the lower boundary of the flow
domain up to 2.3 m or down to 2.7 m below the soil surface (Fig. 2).
Simulation results showed that this did not affect the arrival timeof the tracers. The front and back (receding) parts of the BTCs
were slightly steeper, while the maximum concentration in the well
was about 10% larger for the thicker fi fth layer (bottom boundary
at 2.7 m) than with a 2.3-m-deep bottom boundary.

Water and Solute Fluxes at the Well Boundary
Using HYDRUS-2D we were able to estimate the different waterand solute fluxes through the soil–well boundary. Good agreement(R 2 = 0.944, MIA = 0.961) was obtained between the cumulative water volume flowing to the well (the sum of the cumulative
seepage flux and the variable flux at the bottom of the well) andthe cumulative pumping volume (Fig. 6) during the experiment.
The flux through the seepage face in these simulations started ataround 1.2 d after initiating irrigation. The cumulative seepage
volume first exceeded the cumulative variable flux (groundwater)
after 4.7 d, and at the end of the experiment was about 20% larger
than the variable flux. The time period from 1.5 to 3 d, when mostof the increase in the seepage flux occurred, corresponds roughlywith the time period when the water content at depths below 0.8
m increased as registered with the MCPs (Fig. 4).
The estimated length of the seepage face at the time of the peakconcentrations was about 35 cm. Thus, if a seepage zone exists, it
should have a substantial impact on the solute mass flux (Fig. 7).
Good agreement was obtained between the calculated cumulative
fluxes and the mass recoveries of Br (R 2 = 0.930, MIA = 0.964)
and PFBA (R 2 = 0.950, MIA = 0.969), while the actual mass
recovery of Cl significantly exceeded the simulated value after 4

Fig. 6. Observed and calculated cumulative water fluxes to the well.

d. This may have been caused by the presence of Cl sources in theunsaturated zone that were not accounted for in the model (seeGish and Kung, 2007). All tracers showed much larger simulated
mass fluxes through the estimated seepage face than those throughthe variable boundary at the bottom of the well. Hence, solutes
appeared to be moving mostly along flow paths close to the watertable, then entering the well from the saturated region near thewater table. The seepage face is part of the saturated zone (Fig. 8a),where most of the flow is horizontal (Fig. 8b), and is also the firstpart of the saturated zone encountered by the tracers (Fig. 8c).
Figure 8 shows that the flow patterns in the capillary fringe justabove the phreatic surface are also predominantly horizontal, thus
shortening tracer transport to the well boundary. These findings
are consistent with both recent laboratory (Silliman et al., 2002)
and field observations (Abit et al., 2007). Gardner (2005), in modeling the dynamics of pore water seepage from marsh sediments,
estimated that up to one-third of the total creek-bank discharge
occurred through the seepage face that developed above the tideline at tidal stages below mean tide. Recall that in the submerged
part of the well boundary we used the variable flux boundary condition instead of the equilibrium pressure head distribution. Witha specified pressure head in the submerged zone, the maximum
velocity will occur at the water line in the well and then die off
with depth in the submerged zone (Wilson and Gardner, 2006),
rather than in a more or less uniform fashion as shown in Fig. 8b
for a specified uniform flux. We suggest that such a velocity field
will tighten the concentration contours shown on Fig. 8c, thusincreasing solute flux through the seepage face.
Calculated first arrival times of the tracers to the well were around40 and 30% earlier (for Br and PFBA, respectively) through thevariable flux boundary than through the active seepage face, whichbecame active 1.2 d after irrigation started (Fig. 6). Calculated

www.VadoseZoneJournal.org | 693

Fig. 7. Observed and calculated cumulative solute fluxes of (a) Br, (b) pentafluorobenzoic acid (PFBA), and (c) Cl to the well.

Fig. 8. Simulated distributions of (a) pressure head, (b) fluid flow velocity vectors (actual vector lengths near the well were decreased to fit the figure
scale), and (c) Br concentrations 3 d after irrigation started; SF refers to the seepage face.

cumulative mass fluxes through the seepage face at the end of thesimulation period were 8, 4, and 11 times those through the variable flux boundary (the submerged portion of the well screen) for
Br, Cl, and PFBA, respectively. The highest ratio between thesefluxes was for PFBA because this tracer was applied 7 h after irrigation was initiated and the Cl was applied. As a result, PFBA arrived
at the water table later than the other tracers (i.e., at times when
the seepage face was developing and the seepage flux increased).

www.VadoseZoneJournal.org | 694

The smallest ratio between the seepage solute flux and the variable
flux was for Cl due to its immediate input from the saturated zone
having a background concentration of 7.3 mg L−1. Our results arevery much in agreement with sand box experimental and modeling results by Li et al. (2007), who studied the impact of a seepage
face at an unconfined aquifer–lake interface. They found that mostof the groundwater flow and pollutant flux discharged through anarrow portion near the top of the seepage face.

Well Capture Zone
To avoid solute losses with lateral flow away from the well, theradius of a tracer surface application should be smaller than theradius of the well capture zone. The results of simulations with
the HYDRUS-2D code indicate that the water table divide was
located at about 620 cm from the well. An approximate estimateof the well capture radius can also be obtained with Eq. [11], whichassumes steady-state flow conditions. Such conditions probably
occurred during Days 3 to 5 of the experiment when the pumping rate varied very little around a mean value of 0.12 L s−1. The
assumption of steady-state flow was further confi rmed by the
HYDRUS-2D simulations. Using Eq. [11] and [12], we estimatedthe well capture radius (r*) of our field tracer experiment to be
580 cm and the pumping discharge rate (Qr0) to be 0.12 L s−1.
For the calculations, we used r 0 = 2.5 cm, R = 3000 cm, a = 2000
cm, K = 52 cm h−1, i = 0.41 cm h−1, h0 = 90 cm, and hR = 82 cm.
The values of both r* and Qr0 decreased with increasing hydraulicconductivity (Fig. 9).
The analytical solutions for r* (Eq. [11]) and Qr0 (Eq. [12]) wereobtained with several simplifying assumptions such as steadyflow, having a homogeneous aquifer, and no seepage face. Th iscaused Eq. [12] to underestimate the well capture radius compared with the more comprehensive HYDRUS-2D hydrodynamic
model. Nevertheless, this example demonstrates that for practical
purposes the analytical solution can provide a useful first approximation of system parameters needed for the design of these typesof tracer experiments.

6 Conclusions
The HYDRUS-2D simulations of the field tracer experiment (Gishand Kung, 2007), in which three surface-applied tracers moved
through the unsaturated zone to a pumping well in a perched shallow aquifer, indicate the presence of a seepage face boundary in thewell near the water table. Transport simulations were performed
for several scenarios with different aL values. Values of aL varied
from 5 to 40 cm, while keeping the aL/aT ratio at 10. The best overall results were obtained with a longitudinal dispersivity of 15 cm.
The observed tracer breakthrough curves exhibited steeper risingfronts, more pronounced asymmetry, earlier concentration peaks,
and longer tails than the simulated curves without a seepage face.
Including a seepage face improved the predictions by causing earlier
breakthrough times and higher solute concentrations in the well.

Fig. 9. Well capture radius (r*) and pumping rate (Qr0) as a function of
hydraulic conductivity (Ks) for a given set of parameters (radius of thewell [r0], 2.5 cm; radius of the outer boundary [R], 3000 cm; radiusof the infiltration area, 2000 cm; recharge rate, 0.41 cm h−1; watertable height at radial coordinate r = r0, 90 cm; water table height at r
= R, 82 cm).

Good agreement was obtained between the observed and simulated cumulative water and solute fluxes to the well. The simulation
results indicated that the cumulative water flux to the well throughthe seepage face at the end of the experiment was about 20% larger
than the cumulative variable flux through the saturated zone.
The calculated cumulative seepage face fluxes of the tracers wereapproximately 8, 4, and 11 times the variable flux boundary for
Br, Cl, and PFBA, respectively. These significant differences canbe partially explained by the presence of higher solute concentrations in the capillary fringe and in the upper part of the saturated
zone compared with the lower saturated zone (connected to thesubmerged portion of the well). The higher concentrations fromthe soil solution discharged to the well through the seepage face of
the unconfined aquifer. Our study further suggests that a seepage
face should be considered when evaluating solute fluxes using thismass flux well transport approach.
It must also be realized that an alternative conceptual model, e.g.,simulating preferential flow through a high-hydraulic-conductivity
layer below the water table, can provide a similar or an even better fit
to the observed BTC. Using the preferential flow and deterministictransport model, however, requires knowledge of additional parameters that cannot be assessed with the existing data set. We believethat such a model would also predict the development of the seepage
face; however, the ratio between the solute fluxes through the seepage face and through submerged part of the well would be different
from that obtained with the conventional model.

Acknowledgments

This work has been partially supported through the Interagency USDA and U.S.
NRC Agreement RES-02-008 “Model Abstraction Techniques for Soil Water
Flow and Transport.”

www.VadoseZoneJournal.org | 695

References
Abit, S.M., A. Amoozegar, M.J. Vepraskas, and C.P. Niewoehner. 2007. Solutetransport in the capillary fringe and shallow groundwater: Field evaluaon. Vadose Zone J. 7:890–898.
Anderson, M.P. 1979. Using models to simulate the movement of contaminants through groundwater flow systems. Crit. Rev. Environ. Control
9:97–156.
Aravin, V.I., and S.N. Numerov. 1965. Theory of fluid flow in undeformableporous media. Israel Prog. for Sci. Transl., Jerusalem.
Bear, J. 1972. Dynamic of fluids in porous media. Elsevier, New York.
Buczko, U., and H.H. Gerke. 2006. Modeling two-dimensional water flow andbromide transport in a heterogeneous ligni c mine soil. Vadose Zone J.
5:14–26.
Chenaf, D., and R.P. Chapuis. 2007. Seepage face height, water table posi on,
and well eﬃciency at steady state. Ground Water 45:168–177.
Domenico, P.A., and F.W. Schwartz. 1998. Physical and chemical hydrogeology. 2nd ed. John Wiley & Sons, New York.
Gardner, L.R. 2005. The role of geomorphic and hydraulic parameters ingoverning pore water seepage from salt marsh sediments. Water Resour.
Res. 41:W07010, doi:10.1029/2004WR003671.
Gelhar, L.W., A. Mantaglou, C. Welty, and K.R. Rehfeldt. 1992. A cri cal review
of field-scale dispersion in aquifers. Water Resour. Res. 28:1955–1974.
Gerke, H.H., J. Dusek, T. Vogel, and J.M. Köhne. 2007. Two-dimensional dualpermeability analyses of a bromide tracer experiment on a le-drained
field. Vadose Zone J. 6:651–667.
Gish, T.J., and K.-J.S. Kung. 2007. Procedure for quan fying a solute flux to ashallow perched water table. Geoderma 138:57–64.
Gish, T.J., K.-J.S. Kung, D.C. Perry, J. Posner, G. Bubenzer, C.S. Helling, E.J. Kladivko, and T.S. Steenhuis. 2004. Impact of preferen al flow at varying irriga on rates by quan fying mass fluxes. J. Environ. Qual. 33:1033–1040.
Hodne , M.G., and J. Tomasella. 2002. Marked diﬀerences between van
Genuchten soil water-reten on parameters for temperate and tropicalsoils: A new water-reten on pedo-transfer func ons developed for tropical soils. Geoderma 108:155–180.
Kladivko, E.J., L.C. Brown, and J.L. Baker. 2001. Pes cide transport to subsurface le drains in humid regions of North America. Crit. Rev. Environ. Sci.
Technol. 31:1–61.
Köhne, J.M., and H.H. Gerke. 2005. Spa al and temporal dynamics of preferen al tracer movement towards a le drain. Vadose Zone J. 4:79–88.
Köhne, J.M., B.P. Mohanty, and J. Šimůnek. 2005. Inverse dual-permeability
modeling of preferen al water flow in a soil column and implica ons forfield-scale solute transport. Vadose Zone J. 5:59–76.
Köhne, S., B. Lennartz, J.M. Köhne, and J. Šimůnek. 2006. Bromide transport at a le-drained field site: Experiment, and one- and two-dimensional equilibrium and non-equilibrium numerical modeling. J. Hydrol.
321:390–408.
Kung, K.-J.S., E.J. Kladivko, T.J. Gish, T.S. Steenhuis, G. Bubenzer, and C.S.
Helling. 2000. Quan fying preferen al flow by breakthrough of sequenally applied tracers: Silt loam soil. Soil Sci. Soc. Am. J. 64:1296–1304.
Larsson, M.H., N.J. Jarvis, G. Torstensson, and R. Kasteel. 1999. Quan fying the impact of preferen al flow on solute transport to le drains in asandy field soil. J. Hydrol. 215:116–134.
Legates, D.R., and G.J. McCabe, Jr. 1999. Evalua ng the use of “goodness-offit” measures in hydrologic and hydroclima c model valida on. Water
Resour. Res. 35:233–241.
Li, Y., C. Wang, L.-Z. Yang, and Y.-P. Li. 2007. Influence of seepage face obliquity on discharge of groundwater and its pollutant into lake from a typical unconfined aquifer. J. Hydrodyn. Ser. B. 19:756–761.
Maraqa, M.A., R.B. Wallace, and T.C. Voice. 1997. Eﬀects of degree of water
satura on on dispersivity and immobile water in sandy soil columns. J.
Contam. Hydrol. 25:199–218.
Mohanty, B.P., J.M. Bowman, J.M.H. Hendrickx, J. Šimůnek and M.Th. van
Genuchten. 1998. Preferen al transport of nitrate to a le drain in anintermi ent-flood-irrigated field: Model development and experimental
evalua on. Water Resour. Res. 34:1061–1076.
Mualem, Y. 1976. A new model predic ng the hydraulic conduc vity of unsaturated porous media. Water Resour. Res. 12:513–522.
Neuman, S.P. 1975. Galerkin approach to saturated–unsaturated flow in porous media. p. 201–217. In R.H. Gallagher et al. (ed.) Finite elements in
fluids. Vol. 1. Viscous flow and hydrodynamics. John Wiley & Sons, London.
Nützmann, G., S. Maciejewski, and K. Joswig. 2002. Es ma on of water satura on dependence of dispersion in unsaturated porous media: Experiments and modelling analysis. Adv. Water Resour. 25:565–576.
Petersen, C.T., J. Holm, C.B. Koch, H.E. Jensen, and S. Hansen. 2002. Movement of pendimethalin, ioxinyl and soil par cles to field drainage les.
Pest Manag. Sci. 59:85–96.
Rushton, K.R. 2006. Significance of a seepage face on flows to wells in unconfined aquifers. Q. J. Eng. Geol. Hydrogeol. 39:323–331.

Schuh, W.M., and R.L. Cline. 1990. Eﬀect of soil proper es on unsaturated
hydraulic conduc vity pore-interac on factors. Soil Sci. Soc. Am. J.
54:1509–1519.
Shamsai, A., and T.N. Narasimhan. 1991. A numerical inves ga on of freesurface–seepage face rela onship under steady state flow condi ons.
Water Resour. Res. 27:409–421.
Silliman, S.E., B. Berkowitz, J. Šimůnek, and M.Th. van Genuchten. 2002.
Fluid flow and solute migra on within the capillary fringe. Ground Water
40:76–84.
Simpson, M.J., T.P. Clement, and T.A. Gallop. 2003. Laboratory and numerical inves ga on of flow and transport near a seepage-face boundary.
Ground Water 41:690–700.
Šimůnek, J., and J.A. de Vos. 1999. Inverse op miza on, calibra on and valida on of simula on models at the field scale. p. 431–445. In J. Feyen and
K. Wiyo (ed.) Modelling transport processes in soils at various scales inme and space. Wageningen Pers, Wageningen, the Netherlands.
Šimůnek, J., M. Sejna, and M.Th. van Genuchten. 1999. The HYDRUS-2D
so ware package for simula ng the two-dimensional movement of water, heat, and mul ple solutes in variably-saturated media. Version 2.0.
IGWMC-TPS-53. Int. Ground Water Modeling Ctr., Colorado School ofMines, Golden.
Šimůnek, J., M.Th. van Genuchten, and M. Sejna. 2005. The HYDRUS-1D so ware package for simula ng the one-dimensional movement of water,
heat, and mul ple solutes in variably-saturated media. Version 3.0. Dep.
of Environ. Sci., Univ. of California, Riverside.
Vachaud, G., and M. Vauclin. 1975. Comments on “A numerical model basedon coupled one-dimensional Richards and Boussinesq equa ons.” Water
Resour. Res. 11:506–509.
van Genuchten, M.Th. 1980. A closed-form equa on for predic ng the hydraulic conduc vity of unsaturated soils. Soil Sci. Soc. Am. J. 44:892–898.
Villholth, K.G., and K.H. Jensen. 1998. Flow and transport processes in a macroporous subsurface-drained glacial ll soil: II. Model analysis. J. Hydrol.
207:121–135.
Vogel, T., H.H. Gerke, R. Zhang, and M.Th. van Genuchten. 2000. Modelingflow and transport in a two-dimensional dual-permeability system with
spa ally variable hydraulic proper es. J. Hydrol. 238:78–89.
Wilson, A.M., and L.R. Gardner. 2006. Tidally driven groundwater flow andsolute exchange in a marsh: Numerical simula ons. Water Resour. Res.
42:W01405, doi:10.1029/2005WR004302.

www.VadoseZoneJournal.org | 696