- Soil Science Society of America
To develop predictive models for solute transport in the subsurface, the mechanisms governing transport must be understood and quantified. 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 simulations of the data indicated the presence of an active seepage face along the soil–well interface near the water table. The calculated cumulative water flux to the well through the seepage face was found to be 120% of the variable flux boundary into the well through the submerged zone. In addition, calculated mass fluxes of Br, Cl, and PFBA through the seepage interface were approximately 8, 4, and 11 times, respectively, those through the variable flux boundary. Calculations suggest that a seepage face may be responsible for causing the early arrival of the solutes and the overall shape of the BTCs. Our study indicates the potentially important role a seepage face along the soil–well interface may have in modeling water flow and chemical transport to wells in variably saturated, unconfined aquifers.
Simulation of groundwater flow and solute migration are important for understanding 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 for assessing field-scale water and chemical fluxes through the unsaturated and saturated zones (Gish and Kung, 2007; Gerke et al., 2007).
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 the unsaturated 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 atmospheric pressure 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, a lake, or a pumping well). While seepage faces are often ignored in simulations of flow and transport, they are believed to be important when the thickness of the unsaturated zone is comparable to that of the unconfined 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 of radial flow and transport, Simpson et al. (2003) found that the flow velocity through the seepage face was about 30% higher than at the base of the aquifer, while the solute transport patterns were strongly influenced by the flow characteristic of the system. They concluded that numerical models ignoring the presence of a seepage face may incorrectly predict solute residence time distribution.
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 solute concentration BTCs were collected from a tile drain. This tile drain protocol approach allows chemical transport processes to be precisely modeled during both unsaturated and saturated flow conditions (Kladivko et al., 2001; Petersen et al., 2002; Gish et al., 2004; 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 to a 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 surface-applied tracers. Solute transit times based on a piston flow model predicted that the tracers would reach the water 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, another explanation 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 results with 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 infiltrating 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) to quantify 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).
Materials and Methods
Gish and Kung (2007) presented a detailed description of their field 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 field that 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.
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 was inserted 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 was inserted 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 water contents 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 around the 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 as a 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 end of 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 the observed 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 the solute transport model.
The tracer experimental data were analyzed using the HYDRUS-2D model simulating variably saturated water flow and solute 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 and decay), these equations are, respectively:where θ is the volumetric water content [L3 L−3], t is time [T], ψ is the pressure head [L], K is the hydraulic conductivity tensor [L T−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.  were described using the van Genuchten (1980) relationships:where Se is relative saturation, θr is the residual water content [L3 L−3], θs is the water content at saturation [L3 L−3], Ks is the saturated hydraulic conductivity [L T−1], α [L−1] and n are shape parameters, τ = 0.5, and m = 1 − 1/n. The density effect was not accounted for because measured concentrations of solutes did not exceed 1 g L−1.
Application of Eq. [1–4] to the field tracer experiment using 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 axisymmetrical flow 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 and solute transport. The grid origin (radius r = 0, depth z = 0) was set at the center of the observation well and at a depth of 250 cm from ground level. The finite element mesh size was 5 cm in the vertical direction and 10 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 both flow and solute transport. Along the right boundary, 30 m from the 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 with concentration 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 for radii >2000 cm. Flux concentrations between 10 and 300 cm near the well where the tracers were applied were calculated such that they 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 infiltrating water of 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 (Gish and 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 the well. Although the model stipulated a potential seepage boundary from the elevation of 85 cm to the soil surface, the “active” seepage face occurs only near the water table (see Fig. 2), depending on the transient flow conditions during the experiment. During simulation, the saturated part (ψ ≥ 0) of a potential seepage face is 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 the 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 on transient 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 transient pressure head values. Another reason why we used the variable flux boundary instead of a prescribed pressure head is that we wanted to preserve water fluxes pumped from the well. Note that we did not need 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 the seepage flow and the variable flux), which significantly exceeded the total volume pumped. We gradually decreased the imposed values 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 profile 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 the saturated zone before irrigation.
Calculated Tracer Concentrations 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:where rw is the well radius [L], hw is the time-averaged mean water level 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 the seepage face and variable flux boundaries, respectively, as calculated with HYDRUS-2D, and Qp 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 well and the large water and solute fluxes.
The goodness-of-fit of the model was evaluated using the coefficient of determination (R2) and the modified index of agreement (MIA) as given by Legates and McCabe (1999):where Cw(ti) and Cwobs(ti) represent the simulated and observed concentrations of a tracer in pumped water at time ti, and is the 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 Solution 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 software packages, respectively, as well as laboratory measurements of the saturated water content (θs) of the various soil horizons. During a first step, we used HYDRUS-1D to estimate the hydraulic parameters from the observed water content values that were monitored at different depths in the soil profile 11 m away from the pumping well. By using HYDRUS-1D we assumed 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 compared with 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 water presumably flowed both vertically and horizontally.
For the HYDRUS-1D simulations, we assumed a 250-cm-deep soil profile 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 the fact that markedly different water contents were observed during steady infiltration 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 the flow model, which in total contained 25 hydraulic parameters (the five parameters θr, θs, α, n, and Ks for each of the five layers). To minimize issues of uniqueness in the inverse solution, we decreased the number of optimized parameters by using experimentally measured 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 three parameters (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 was low, and no rain or irrigation occurred in the week before the experiment (Gish and Kung, 2007). Therefore, we assumed that the initial pressure head profile was close to equilibrium with the water 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 around the well after irrigation began. As a result, some lateral water flow in the saturated zone occurred, which had to be accounted for by the flux through the lower boundary in the one-dimensional HYDRUS-1D simulations. This flux varied with time, depending on the position of the water table. To simulate this process, we used the HYDRUS-1D horizontal drain boundary condition with a drain 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 most closely represented our field conditions.
Initial estimates of the α and n parameters for the inverse procedure were based on the assumed initial equilibrium distribution of the pressure head and the observed water contents before the experiment, while Ks values were initially estimated from soil textural class and particle size distribution information. We used a sequential inverse procedure by starting the parameter search with the first layer, while assuming that the other parameters were known (the initial estimates) and using observed water contents of 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.3 m, 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 procedure by gradually increasing the number of layers and hence the number of 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 Conditions
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 table divide, where the radial flow velocity is equal to zero (Fig. 3 ). The radius of the water table divide can be roughly estimated from a solution of the Boussinesq equation for steady-state radial unconfined flow toward a well, subject to Dirichlet boundary conditions as follows: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 infiltration area (r0 ≤ a ≤ R), r0 and R are the radii of the well and the outer boundary [L], respectively, h0 and hR are the water table heights at r = r0 and r = R, respectively, and r is the radial coordinate.
Assuming a constant hydraulic conductivity, Ks, and integrating Eq. , we obtained an analytical expression for the radial water flux to the well, Q:For the limiting case of a = R, this solution coincides with an equation derived by Aravin and Numerov (1965).
An equation for the water table divide radius (r*) may be obtained by solving the equation for Q = 0 for the case when r* < a:The pumping well discharge is calculated from
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 R2 = 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 profile.
Table 1 presents the fitted hydraulic parameters obtained using the invoked inverse procedure. The values of α and n are within the range reported by Hodnett and Tomasella (2002) for soils with similar textures. A relatively high value of 10 was obtained for the van Genuchten hydraulic parameter n for the fifth layer (gravely sand), 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. While some 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 result of inverse simulations for the limited range of water content in the unsaturated zone and for the most common value of the pore-connectivity parameter τ = 0.5 (Eq. ). The value of τ can vary in a range from −1 to 2.5 (Mualem, 1976), however, and Schuh and Cline (1990) reported that τ varied between −8.73 and 14.80. This most significantly affects the hydraulic conductivity values at water 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 two-dimensional axisymmetrical flow model against the observed water contents. 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 (gravely sand), 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 an imperfect 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 R2 = 0.858, slightly smaller than the value of 0.877 obtained with the HYDRUS-1D calibration.
Despite the achieved good fit to the experimental data, we realize that 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 and then 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 unique soil hydraulic parameters. This is reasonable since no attempts were made to extrapolate the data observed beyond the experimental time period.
Two-Dimensional Simulations of Water Flow and 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 the field tracer tests (Gelhar et al., 1992). We used the same value for all layers, while assuming in most cases a ratio of 10 between the longitudinal (aL) and transversal (aT) dispersivities. This ratio is well within the range of ratios (5–20) reported in several studies by Anderson (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 the dispersivity has not yet been established, however, and therefore we used constant values of these parameters for the entire variably saturated domain. Figure 5 compares the observed and calculated (using Eq. ) solute BTCs of the tracers in the well. The results are shown for 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 BTCs of 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. ) 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 much earlier arrival times and lower maximum concentrations than the observed values. Decreasing the ratio aL/aT from 10 to 5 did not have 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). 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 (the pumping rate generally increased with time) and because the pump had to be replaced at some point during the experiment (Gish and Kung, 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 these fluxes 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 found to produce similar results. We additionally tested the effect of having a different thickness of the fifth layer (where most of the well 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 time of 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 fifth 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 water and solute fluxes through the soil–well boundary. Good agreement (R2 = 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) and the cumulative pumping volume (Fig. 6 ) during the experiment. The flux through the seepage face in these simulations started at around 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 most of the increase in the seepage flux occurred, corresponds roughly with 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 peak concentrations 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 (R2 = 0.930, MIA = 0.964) and PFBA (R2 = 0.950, MIA = 0.969), while the actual mass recovery of Cl significantly exceeded the simulated value after 4 d. This may have been caused by the presence of Cl sources in the unsaturated zone that were not accounted for in the model (see Gish and Kung, 2007). All tracers showed much larger simulated mass fluxes through the estimated seepage face than those through the variable boundary at the bottom of the well. Hence, solutes appeared to be moving mostly along flow paths close to the water table, then entering the well from the saturated region near the water 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 first part of the saturated zone encountered by the tracers (Fig. 8c). Figure 8 shows that the flow patterns in the capillary fringe just above 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 tide line 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. With a 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, thus increasing solute flux through the seepage face.
Calculated first arrival times of the tracers to the well were around 40 and 30% earlier (for Br and PFBA, respectively) through the variable flux boundary than through the active seepage face, which became active 1.2 d after irrigation started (Fig. 6). Calculated cumulative mass fluxes through the seepage face at the end of the simulation 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 these fluxes 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).
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 are very 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 most of the groundwater flow and pollutant flux discharged through a narrow portion near the top of the seepage face.
Well Capture Zone
To avoid solute losses with lateral flow away from the well, the radius of a tracer surface application should be smaller than the radius 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 estimate of the well capture radius can also be obtained with Eq. , which assumes 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 confirmed by the HYDRUS-2D simulations. Using Eq.  and , we estimated the 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 r0 = 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 hydraulic conductivity (Fig. 9 ).
The analytical solutions for r* (Eq. ) and Qr0 (Eq. ) were obtained with several simplifying assumptions such as steady flow, having a homogeneous aquifer, and no seepage face. This caused Eq.  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 types of tracer experiments.
The HYDRUS-2D simulations of the field tracer experiment (Gish and 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 the well 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 rising fronts, 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.
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 through the 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 were approximately 8, 4, and 11 times the variable flux boundary for Br, Cl, and PFBA, respectively. These significant differences can be 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 the submerged portion of the well). The higher concentrations from the 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 this mass 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 deterministic transport model, however, requires knowledge of additional parameters that cannot be assessed with the existing data set. We believe that 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.
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.”
All rights reserved. No part of this periodical may be reproduced or transmitted in any form or by any means, electronic or mechanical, including photocopying, recording, or any information storage and retrieval system, without permission in writing from the publisher.
- Received May 5, 2009.