Advances in Water Resources 26 (2003) 989–1000 www.elsevier.com/locate/advwatres
Unsaturated hyporheic zone flow in stream/aquifer conjunctive systems Garey A. Fox
a,*
, Deanna S. Durnford
b
a b
Department of Civil Engineering, University of Mississippi, University, MS 38677-1848, USA Department of Civil Engineering, Colorado State University, Fort Collins, CO 80523-1372, USA
Received 15 February 2002; received in revised form 13 February 2003; accepted 15 May 2003
Abstract Saturated flow is typically assumed for seepage from a stream underlain by an alluvial aquifer. However, if the water table falls a sufficient distance below a semipervious streambed, the head losses in this less conductive layer will cause the region beneath the stream, or hyporheic zone, to become unsaturated. Hyporheic zone flow is defined loosely in this research as the flow that occurs underneath the streambed. Unsaturated flow transforms streams from constant head boundaries to constant flux boundaries, impacting the biogeochemistry in the hyporheic zone. The objective of this paper is to discuss the development and implications of unsaturated flow beneath the streambed. Conditions under which saturated or unsaturated flow occurs and the characteristics of each flow regime are discussed. Next, the effect of unsaturated flow is illustrated for the case of stream leakage induced by a well pumping from an aquifer that is hydraulically interacting with a partially penetrating stream. Prior analytical solutions for alluvial well depletions fail to model unsaturated flow between the streambed and water table. An approximating solution is proposed to estimate aquifer drawdown and stream depletion under saturated/unsaturated hyporheic zone flow conditions. 2003 Elsevier Ltd. All rights reserved. Keywords: Stream/aquifer interaction; Hyporheic zone; Unsaturated flow; Stream leakage; Constant head boundary; Constant flux boundary
1. Introduction The hydraulic connection between surface water and groundwater resources has long been recognized as a crucial factor in the management of water resources. However, new questions in water supply, water quality, and stream ecology are emerging as water supplies become overappropriated, water quality degradation becomes more widespread, and water management must satisfy the dual objectives of water supply and ecosystem habitat conservation [1,29]. These new concerns are causing researchers and managers to take a new look at the complex hydrology at the interface between surface water and groundwater. An important groundwater/surface water interface is stream leakage to an underlying aquifer. Unsaturated flow can occur within this interface under a number of different scenarios. Ephemeral streams are commonly perched above the groundwater table. In cases where the *
Corresponding author. E-mail addresses:
[email protected] [email protected] (D.S. Durnford).
(G.A.
Fox),
durn-
0309-1708/$ - see front matter 2003 Elsevier Ltd. All rights reserved. doi:10.1016/S0309-1708(03)00087-3
stream and aquifer are initially connected, groundwater extraction from wells can initiate unsaturated flow. In this case, initiation of unsaturated flow is dependent on the streambed hydrologic properties. If the streambed is thin with approximately the same physical and hydraulic properties as the underlying soil, flow through the hyporheic zone is saturated, and the discharge will be a function of the location of the water table. However, the streambed often has a conductivity that is orders of magnitude less than the underlying aquifer and is sufficiently thick so that most of the hydraulic head loss occurs in the less conductive layer [7,21,23,25]. This conductivity contrast is more realistically analyzed as a layered system rather than a homogeneous profile. In general, seepage from a stream with a less permeable streambed can be characterized by three hyporheic zone flow regimes: saturated flow, a transition zone, and unsaturated flow. The location of the water table, soil properties, and accessibility of an air phase determine the regime. The objective of this paper is to discuss the development and implications of unsaturated hyporheic zone flow. First, the conditions that determine the flow
990
G.A. Fox, D.S. Durnford / Advances in Water Resources 26 (2003) 989–1000
Fig. 1. Outline of the problem considered by (a) Theis [26], (b) Hantush [18], and (c) Hunt [19]. L ¼ distance between the stream and pumping well and Q ¼ discharge rate of the pumping well. Modified from [19].
regime and the characteristics of each regime are discussed. Then, the effect of an unsaturated hyporheic zone is illustrated for the case of stream leakage induced by a well pumping from an aquifer that is hydraulically interacting with a partially penetrating stream. A pumping well will intercept groundwater, reducing recharge to the stream and/or extracting water from the stream. While this problem has been analyzed extensively for decades [6,9,18–20,25,26], it is always assumed that the stream is a constant head boundary and flow between the stream and aquifer is saturated. The most well known equation for drawdown due to a pumping well in an infinite aquifer is the Theis solution [27]. Theis expanded this equation to simulate a pumping well adjacent to a stream through image well theory [26]. This solution assumes an infinitely long, straight, completely penetrating stream in a homogeneous, isotropic, confined aquifer stressed by a well pumping at a constant discharge rate, as shown in Fig. 1a. Glover and Balmer [17] expressed the Theis solution [26] in terms of stream depletion. Hantush [18] developed an analytical model that considers the effects of a semipervious streambed. The semipervious streambed is represented as a vertical layer of lower conducting material extending throughout the saturated thickness of the aquifer (Fig. 1b). Hunt [19] developed an analytical solution for drawdown and stream depletion that accounts for a semipervious streambed and stream partial penetration (Fig. 1c). Butler et al. [6] developed an analytical solution that considers finite stream width and an aquifer of limited lateral extent. However, these analytical solutions fail to model the case where the region between the streambed and water table becomes unsat-
urated. If the water table falls a sufficient distance below the bottom of the streambed, the higher head loss in the less conductive streambed will cause the hyporheic zone to become unsaturated. When this occurs, the stream and aquifer are often said to hydraulically ‘‘disconnect’’ [4,24]. This term is a misnomer. The stream and underlying aquifer remain hydraulically connected, albeit through unsaturated flow. A more appropriate term is that the stream is ‘‘perched’’ above the water table. Unsaturated flow is important mathematically because it transforms streams from constant head boundaries, where the discharge from the stream is dependent on the drawdown in the aquifer, to constant flux boundaries, where the discharge is independent of the position of the water table. Unsaturated hyporheic zone flow is also important environmentally because it impacts biogeochemical reactions in the hyporheic zone.
2. Stream/aquifer interaction A homogeneous and isotropic streambed is assumed to exist between the stream and underlying hyporheic zone. The streambed is assumed to have a lower hydraulic conductivity than the underlying soil. It is also assumed that leakage from the stream to the aquifer is steady, vertical, one-dimensional flow. For these conditions, it can be shown that the streambed will remain saturated [10]. The following discussion addresses three stream/aquifer hydrologic states or regimes: • Regime A––Saturated flow: Stream leakage rate depends on the location of the water table.
G.A. Fox, D.S. Durnford / Advances in Water Resources 26 (2003) 989–1000
• Regime B––Transition zone: Stream leakage rate governed by unsaturated flow hydraulic conditions. • Regime C––Unsaturated gravity-driven flow: Stream leakage rate is not a function of the location of the water table. 2.1. Regime A: saturated flow Saturated flow occurs in the hyporheic zone when the water table is located within the streambed, when the water table is slightly below the bottom of the streambed but water pressures are not sufficiently negative to desaturate the subsoil, or if there are no pathways for the air phase to enter the hyporheic zone. When the water table is located at an elevation within the streambed, the specific discharge, q, through that region of the streambed between the bottom of the stream and the water table (assuming one-dimensional, vertical flow, where a downward flux is assumed to be negative) is sw q ¼ Ksb ð1Þ sw Hw where Ksb is the saturated hydraulic conductivity of the streambed and sw is drawdown, defined as the distance between the water level in the stream, Hw , and the water table. If the water table is located just at the bottom of the streambed, Eq. (1) reduces to Hw þ M q ¼ Ksb M
ð2Þ
where M is the thickness of the streambed. In both of these simple cases, the hydraulic head in the stream is driving the flow and viscous resistance is proportional to the hydraulic conductivity of the streambed. If the water table declines to an elevation below the bottom of the streambed, the specific discharge is given by q ¼ Ksb
Hw þ M hw M
ð3Þ
where hw is the water pressure head at the bottom of the streambed. Because the water table is below the streambed, the water pressure head, hw , at the streambed/hyporheic zone interface is negative. It is useful to consider water pressure head, hw , in terms of capillary pressure head, hc . If pore fluid pressures are measured with respect to ambient air pressures, assumed zero, then hc ¼ hw
ð4Þ
Rewriting Eq. (3) in terms of capillary pressure head, hc , gives Hw þ M þ hc q ¼ Ksb M
ð5Þ
which is similar to the expression used by Bouwer [4]. The hyporheic zone below the streambed may or may not desaturate, depending on the magnitude of the water
991
pressure head at the interface between the streambed and hyporheic zone. Desaturation of the underlying soil will occur if the capillary pressure head is greater (i.e., the water pressure head is more negative) than an air entry capillary pressure head, he . The air entry pressure head is that capillary pressure head below which the nonwetting phase (in this case air) becomes discontinuous and is not capable of flowing through the material [10]. The entry pressure head loosely represents the height of the capillary fringe in a uniform soil. The value of an entry pressure head can be estimated by curve fitting capillary pressure-saturation data [10,11]. There is considerable ambiguity regarding this critical capillary pressure head and it is variously referred to as a displacement pressure head, a bubbling pressure head, and an entry pressure head [10,11,30]. If there are no pathways for air to enter the hyporheic zone, the region below the streambed will also remain saturated no matter how far the water table drops below the streambed. In this case, water pressure heads could theoretically become very large negative values. This paper will assume that pathways exist for air to enter the hyporheic zone beneath the streambed. 2.2. Regime B: transition flow Regime B is a transition zone between saturated hyporheic zone flow and gravity-driven, unsaturated hyporheic zone flow. Unsaturated hyporheic zone flow occurs when the water table falls a sufficient distance below the bottom of the streambed so that some of the pores desaturate, but the water table is not sufficiently far below the streambed for a unit hydraulic gradient to exist. It can be shown that the streambed will remain saturated [10]. Fig. 2a shows the distribution of capillary pressure head in the hyporheic zone for this regime. The specific discharge of water through the unsaturated hyporheic zone is given by ohc q ¼ Kðhc Þ 1 ð6Þ oz where hc is the capillary pressure head, Kðhc Þ is a constitutive relation between hydraulic conductivity and capillary pressure, and z is the depth within the hyporheic zone. It is convenient to fit algebraic expressions to the capillary pressure head–hydraulic conductivity relationship. Among the most popular of these algebraic expressions are the Brooks–Corey [5,10] and van Genuchten [28] equations. The Brooks–Corey equations are Kðhc Þ ¼ Ks Kðhc Þ ¼ Ks
hc 6 he g he hc > he hc
ð7Þ
where Ks is the saturated hydraulic conductivity, hc is capillary pressure head, he is the entry capillary pressure
992
G.A. Fox, D.S. Durnford / Advances in Water Resources 26 (2003) 989–1000 Table 1 Typical values of entry pressure head, he , and Brooks–Corey parameter, g, on a drainage curve Soil type
he (cm of water)
Coarse sand Medium sand Fine sand Silt Clay
2–10 10–35 35–70 70–150 >200–400
g 8.0 6.5 5.9 3.1 2.6
Values are obtained from Bear [3] and Carsel and Parrish [8].
soils are shown in Table 1. Using these representative values, Eq. (8) is solved for the change in the capillary pressure head gradient with elevation ðohc =ozÞ as a function of the capillary pressure head ðhc Þ. Results are shown in Fig. 3 for fine, medium and coarse subsoils with two different ratios of the specific discharge to saturated hydraulic conductivity, jqj=Ks . When ohc =oz ¼ 0, a unit hydraulic gradient exists in the hyporheic zone. The formation of this unit hydraulic gradient signals the end of regime B. This figure demonstrates that the capillary pressure gradient ohc =oz goes to zero quickly, especially for medium and coarse sand subsoils. Fig. 2. Capillary pressure head ðhc Þ distribution in the hyporheic zone (a) during regime B, prior to the formation of a unit hydraulic gradient in the hyporheic zone and (b) during regime C, after formation of a unit hydraulic gradient. hcL ¼ interface capillary pressure head, Hw ¼ river stage, sw ¼ drawdown, and M ¼ streambed thickness.
head, and g is a parameter dependent on the pore-size distribution index [5,10]. Rearranging Eq. (6) and replacing Kðhc Þ with the Brooks–Corey equation results in g ohc jqj hc ¼1 hc > he ð8Þ Ks he oz
2.3. Regime C: unsaturated gravity-driven hyporheic zone flow Regime C is defined by the presence of only gravitydriven flow ðohc =oz ¼ 0Þ in contrast to the other regimes where there is also a pressure gradient in the unsaturated hyporheic zone beneath the streambed (Fig. 2b). In regime C, stream leakage reaches a maximum as the interface capillary pressure head, hcL , becomes a maximum. This maximum interface capillary pressure
Other constitutive relationships for Kðhc Þ can be used instead of the Brooks–Corey equations. Eq. (8) can be integrated from the water table to the bottom of the streambed to relate the drawdown, sw , and the capillary pressure head at the bottom of the streambed, hcL , in terms of the hydraulic conductivity of the subsoil and soil parameters, he and g Z hcL ohc g sw M Hw ¼ ð9Þ jqj hc 0 1 Ks he Eq. (5) with hc ¼ hcL and Eq. (9) can be solved simultaneously for the specific discharge, q, and the interface capillary pressure head, hcL . The specific discharge in regime B has a magnitude between those in regimes A and C, but under most circumstances, the length of time that this regime applies is small. Representative values of the entry pressure head, he , and Brooks–Corey parameter, g, for different sub-
Fig. 3. Gradient of capillary pressure head with elevation ðohc =ozÞ as a function of the interface capillary pressure head ðhc Þ for different soil types and the ratio of specific discharge to saturated hydraulic conductivity, jqj=Ks .
G.A. Fox, D.S. Durnford / Advances in Water Resources 26 (2003) 989–1000
993
head will be referred to as the ultimate interface capillary pressure head, hcu . This ultimate interface capillary pressure head is a function of the air entry pressure head, he , and saturated hydraulic conductivity, Ks , of the aquifer hcu ¼ he
Ks jqj
1=g ð10Þ
The specific discharge, qmax , is given by ðhcu þ Hw Þ qmax ¼ Ksb 1 þ M
ð11Þ
and simplifying Eq. (6) as qmax ¼ Ks
he hcu
g hcu > he
ð12Þ
Eqs. (11) and (12) can be solved simultaneously for qmax and hcu . In regime C, the specific discharge and ultimate capillary pressure head are independent of the position of the water table. Thus, the stream changes from a constant head boundary to a constant flux boundary. Fig. 4 shows an example assuming a stream with flow depth, Hw ¼ 0:5 m, a 1-m thick streambed, and hydraulic conductivity, Ksb , of 0.1 m/day, overlays a subsoil with a saturated hydraulic conductivity, Ks , of 10 m/day, g ¼ 6:5, and he ¼ 0:1 m. These values are characteristic of a silt streambed overlying a medium, sandy aquifer [8,22]. Fig. 4a plots the specific discharge, nondimensionalized by Ksb , as a function of position of the water table below the streambed, nondimensionalized by the streambed thickness, M. This figure shows that (a) stream leakage varies with water table position in regimes A and B, (b) the transition from regimes A– C occurs over a small range of drawdowns, (c) a constant maximum leakage rate, jqmax j, is obtained when a unit hydraulic gradient occurs in regime C, and (d) the stream leakage does not change with any additional decline in water table position in regime C. Fig. 4b plots the interface capillary pressure head ðhcL Þ, nondimensionalized by the entry pressure head of the subsoil ðhe Þ, as a function of the position of the water table below the streambed. During saturated flow (regime A), pore spaces are completely filled with water because the water pressure heads are not sufficiently negative to desaturate the pores. The interface capillary pressure head increases with increasing drawdown in regime B, reaching a constant, maximum capillary pressure head. This maximum capillary pressure head, or ultimate interface capillary pressure head, hcu , indicates gravity-driven flow is occurring in the unsaturated hyporheic zone. Once hcu is reached, the capillary
Fig. 4. (a) Specific discharge from the stream, jqj, divided by the streambed hydraulic conductivity, Ksb , and (b) interface capillary pressure head, hcL , divided by the entry pressure head of the aquifer, he , as a functions of water table position for regimes A (saturated flow), B (transition), and C (unsaturated flow). Hw ¼ river stage, sw ¼ drawdown, and M ¼ streambed thickness.
pressure head does not increase with additional declines in the water table.
3. Drawdown and stream depletion for saturated/unsaturated hyporheic zone flow Stream depletion and aquifer drawdown in response to a pumping well hydraulically interacting with a stream has been the subject of intensive research for decades. Recently, Hunt [19] developed analytical solutions for drawdown and stream depletion from a partially penetrating stream with a semipervious streambed. However, this solution is limited to saturated flow between the stream and aquifer (i.e., regime A). In the following sections, HuntÕs solution is introduced and then equations are presented that modify his solution to analyze drawdown and stream depletion for both
994
G.A. Fox, D.S. Durnford / Advances in Water Resources 26 (2003) 989–1000
saturated and gravity-driven, unsaturated hyporheic zone flow (regime C). 3.1. Hunt’s [19] solution for saturated flow The governing equation for HuntÕs solution is a twodimensional, partial differential equation with a sink due to pumping, Q, and recharge source from the stream, ksw 2 o s w o2 s w osw Qdðx LÞdðyÞ þ ksw dðxÞ T þ ¼S ox2 oy 2 ot ð13Þ where x and y are horizontal coordinates with respect to a datum at the center of the river on a perpendicular line with the well, t is the time since the start of pumping, T is transmissivity, S is storage coefficient, Q is a constant well discharge rate, d is the Delta Dirac function, L is the distance between the pumping well and the center of the stream, and k is a streambed conductance coefficient given by k¼
Ksb W M
ð14Þ
where Ksb is the streambed hydraulic conductivity, W is stream width, and M is streambed thickness. The third term on the right side of Eq. (13) accounts for stream leakage and explicitly shows a dependency on drawdown. The variable k is derived from the total discharge per unit stream length ðQ0 Þ, which Hunt [19] assumes as s w ð15Þ Q0 ¼ Ksb W ¼ ksw M where sw is the drawdown in the aquifer. The stream is assumed partially penetrating with a semipervious streambed and the aquifer is homogeneous, isotropic and of infinite horizontal extent. The stream is modeled as a constant head boundary with vertical and horizontal streambed cross-sections small compared to the aquiferÕs saturated thickness. Hunt derived an equation for drawdown, sw , in terms of the Theis well function, E1 [27] ( " # Q ðL xÞ2 þ y 2 E1 sw ðx; y; tÞ ¼ 4pT 4Tt=S " # ) Z 1 ðL þ jxj þ 2T h=kÞ2 þ y 2 h e E1 dh 4Tt=S 0
rffiffiffiffiffiffiffi! 2 SL2 kt kL þ exp 4ST 2T 4Tt 0sffiffiffiffiffiffiffiffi rffiffiffiffiffiffiffi1 k2 t SL2 A erfc@ þ 4ST 4Tt
Qs ¼ erfc Q
ð18Þ
HuntÕs model assumes that changes in water table elevation are small, and the stream/aquifer connection is saturated [19,24]. 3.2. Proposed solution for saturated/unsaturated flow When a pumping well is located sufficiently close to a stream and/or the groundwater extraction rate causes significant depletion, part of the stream may become perched above the water table [24]. The stream will initially become perched nearest the well and then the length of stream undergoing unsaturated flow increases over time as pumping continues. Therefore, a more realistic model is one in which segments of the stream closest to the pumping well become perched while other contributing stream segments remain connected to the aquifer by saturated flow, as shown in Fig. 5, where yhd is the time-varying half-length of stream with unsaturated flow. This research proposes an approximating solution to model the scenario presented in Fig. 5, making the following assumptions: (1) changes in the water surface elevation due to pumping are small, (2) vertical and horizontal streambed cross-sections are small compared to the aquiferÕs saturated thickness, (3) ratio of vertical to horizontal flow components in the aquifer is small (Dupuit flow assumption), (4) drawdown is small compared to the saturated thickness of the aquifer, and (5) the transition flow regime can be as-
ð16Þ where E1 ðuÞ ¼
Z u
1
1 x e dx x
ð17Þ
and u is the argument of the well function. Hunt [19] also derived an equation for stream depletion, Qs as a function of the pumping well discharge, Q
Fig. 5. Coordinate system and variable definition for the proposed model. yhd ¼ time-varying length of stream for which unsaturated hyporheic zone flow occurs, L ¼ distance between the stream and pumping well, and Q ¼ discharge rate of the pumping well.
G.A. Fox, D.S. Durnford / Advances in Water Resources 26 (2003) 989–1000
sumed to have a negligible impact on stream/aquifer exchange. The drawdown contributions of flow from the unsaturated flow stream segments and the stream segments that remain saturated are linearly superimposed. Superposition is appropriate since the assumptions result in linear governing partial differential equations. The influence of the groundwater extraction well is simulated using the Theis equation [27] " # 2 Q ðL xÞ þ y 2 E1 sT ¼ ð19Þ 4pT 4Tt=S Hunt [19] gives the equation governing the recharge contribution by an infinite stream undergoing saturated flow as " # Z 1 2 Q ðL þ jxj þ 2T h=kÞ þ y 2 h e E1 dh ð20Þ 4pT 0 4Tt=S where Q is the constant discharge rate of the pumping well. However, during unsaturated flow, segments of the stream closest to the pumping well become perched while other contributing stream segments remain connected through saturated flow. This research hypothesizes that the recharge contribution to drawdown by these saturated flow stream segments can be approximated using a modification of HuntÕs integral expression given by Eq. (20). The modification of HuntÕs integral expression is based on integrating only along the stream reach that remains connected to the aquifer through saturated flow rather than over an infinite stream length. This is accomplished by modifying the lower limit of the integration in Eq. (20), such that " # Z 1 Q ðL þ jxj þ 2T h=kÞ2 þ y 2 h e E1 dh ð21aÞ 4pT h0 4Tt=S
995
with unsaturated flow, and t0 ðlÞ is the time period that unsaturated flow has occurred. This equation is developed assuming unit gradient conditions (i.e., regime C) throughout the unsaturated zone. Eq. (22) is numerically equivalent to GloverÕs analytical solution for the rise of a groundwater mound due to continuous recharge over a rectangular strip [16]. The time period that unsaturated flow has occurred at a particular stream location within yhd depends on l. Unsaturated flow initiates at the stream location nearest the pumping well and then the length of stream with unsaturated flow increases over time. At the stream location nearest to the pumping well, recharge to the alluvial aquifer occurs at the constant rate qmax for the entire time period of unsaturated flow. However, at the stream location dividing saturated and unsaturated flow, recharge is just beginning to occur at the rate qmax . The solution assumes a squared distribution, as shown in Fig. 6, for the relationship between the time period that unsaturated flow has occurred and the location, l, along the stream reach with unsaturated flow, t0 ðlÞ. This distribution is selected based on comparison to an updated MODFLOW RIVER package, RIV_AC [14]. Linear superposition results in an equation for aquifer drawdown when the stream has both saturated and unsaturated flow sw ðx; y; tÞ ¼ sT ss sU
ð23Þ
The first term on the right-hand side of Eq. (23) is equivalent to the Theis [27] equation. The second term is the drawdown contribution due to recharge from stream segments with regions beneath the streambed undergoing saturated flow, as given by Eq. (21a). The last term is the drawdown contribution due to recharge from stream segments with zones undergoing unsaturated flow, given by (22).
where h0 ¼
f ðT ; L; kÞyhd k 2T
ð21bÞ
The function, f ðT ; L; kÞ, in h0 is derived empirically based on comparisons to an updated MODFLOW RIVER package, RIV_AC [14] rffiffiffiffiffiffi T f ðT ; L; kÞ ¼ ð21cÞ Lk The recharge contribution to drawdown from the unsaturated stream segment, sU , is modeling as an infinite line of recharge wells supplying water to the alluvial aquifer at the rate qmax " # Z yhd 2 qmax W x2 þ ðy lÞ E1 sU ¼ dl ð22Þ 4Tt0 ðlÞ=S yhd 4pT where W is the stream width, l is the variable of integration relating to the length along the stream segment
Fig. 6. Approximating relationship for the time period that unsaturated flow occurs, t0 ðlÞ versus location along the stream with unsaturated flow, l.
996
G.A. Fox, D.S. Durnford / Advances in Water Resources 26 (2003) 989–1000
An approximation for the stream depletion occurring during saturated/unsaturated hyporheic zone flow is derived based on the streamflow contribution from the perched stream segment and the streamflow contribution from the saturated flow stream segment. The stream depletion at the moment that the stream begins to perch is 0 sffiffiffiffiffiffiffiffiffiffi! 2 SL2 k thd kL @ Qs ðthd Þ ¼ Q erfc exp þ 4Tthd 4ST 2T 0sffiffiffiffiffiffiffiffiffiffi sffiffiffiffiffiffiffiffiffiffi11 k2 thd SL2 AA ð24Þ erfc@ þ 4ST 4Tthd Once the stream begins to perch (i.e., t > thd ), additional stream depletion is derived from the stream segments with unsaturated hyporheic zone flow. If the recharge wells representing this unsaturated stream segment are assumed to operate for an equivalent time period, the volumetric discharge from this part of the stream ðQUNSAT Þ is given by ðhcu þ Hw Þ QUNSAT ðtÞ ¼ jqmax jWyhd ¼ Ksb 1 þ Wyhd M ð25Þ However, since the unsaturated flow initiates at the stream location nearest the pumping well and then the length of stream with unsaturated flow increases over time, QUNSAT is given by R yhd q Wt0 ðlÞdl yhd max ð26Þ QUNSAT ðtÞ ¼ t thd where t0 ðlÞ is assumed to be governed by a squared distribution as shown in Fig. 6. The stream depletion from saturated flow stream segments ðQSAT Þ is dependent on the drawdown Z 1 QSAT ðtÞ ¼ 2k sw ð0; y 0 ; tÞdy 0 ð27Þ yhd
Note that this equation makes use of the symmetry relationship along the saturated stream segments. The stream depletion can then be expressed as a function of the stream depletion at the onset of unsaturated flow, the stream depletion from the unsaturated stream segment, and the stream depletion from the saturated flow stream segment using numerical techniques Qs ðto Þ ¼ Qs ðthd Þ Qs ðti Þ ¼ Qs ðti1 Þ þ fQUNSAT ðti Þ QUNSAT ðti1 Þg þ fQSAT ðti Þ QSAT ðti1 Þg
ð28Þ
The following steps summarize the solution methodology • Solve (11) and (12) simultaneously for qmax and hcu .
• Solve for the time, thd , at which unsaturated flow initiates at the stream location nearest the well, i.e., ðx; yÞ ¼ ð0; 0Þ, using (16). The time the stream begins to perch is when sw ð0; 0; thd Þ ¼ Hw þ M þ hcu . • If t < thd , then drawdown and stream depletion are given by (16) and (18). • If t > thd , then 1. Solve for yhd iteratively using (23), where sw ð0; yhd ; tÞ ¼ Hw þ M þ hcu . 2. Solve Eqs. (24), (26) and (27) for the components of stream depletion. 3. Evaluate sw ðx; y; tÞ and QS using Eqs. (23) and (28). Since the approximating solution is based on the methodology of Hunt [19], the solution is limited to certain stream/aquifer conditions. The primary limitation is that horizontal flow must predominate in the aquifer. Critical considerations for meeting this criterion include the distance between the pumping well and stream, anisotropy ratio, stream width, and aquifer thickness, as discussed by Butler et al. [6]. A partially penetrating pumping well will create significant vertical flow components at the stream if the distance between the stream and pumping well ðLÞ is less than 1=2 2bðKh =Kv Þ where b is the aquifer thickness, Kv is the vertical aquifer hydraulic conductivity, and Kh is the horizontal aquifer hydraulic conductivity. A moderate degree of anisotropy ðKh =Kv ¼ 10Þ can be neglected if the ratio between the distance between the pumping well and stream ðLÞ to the stream width ðW Þ is greater than 15. For a high degree of anisotropy ðKh =Kv ¼ 100Þ, the impact of anisotropy can be neglected when the ratio of distance between the pumping well and stream ðLÞ to the stream width ðW Þ is greater than 75 [6]. Because the solution superimposes contributions from unsaturated and saturated flow stream segments that are represented by line sources, stream width effects can influence the validity of the approximating solution. Research has indicated that stream width effects become significant when the ratio of the distance between the stream and pumping well ðLÞ to stream width ðW Þ is less than 25–30 [15].
4. Discussion The significance of unsaturated hyporheic zone flow is investigated using a hypothetical stream/aquifer system whose properties are given in Table 2. The aquifer is coarse sand with a hydraulic conductivity 100 times greater than the overlying semipervious streambed layer. Specific yield, Sy , is used in place of the storage coefficient, S, to simulate the storage capabilities of an unconfined aquifer. The proposed saturated/unsaturated approximating solution is compared to the analytical solution of Hunt [19], which assumes a saturated stream/
G.A. Fox, D.S. Durnford / Advances in Water Resources 26 (2003) 989–1000
997
Table 2 Characteristic aquifer, stream, and streambed properties for a hypothetical stream/aquifer system Aquifer parameters Pumping rate, Q (m3 day1 ) Saturated hydraulic conductivity, Ks (m day1 ) Transmissivity, T (m2 day1 ) Specific yield, Sy Entry pressure head, he (m) Brooks–Corey parameter, g Pumping well location Distance between stream and pumping well, L (m) Observation well location x-coordinate location, x (m) y-coordinate location, y (m) Stream and streambed parameters Streambed hydraulic conductivity, Ksb (m day1 ) Stream width, W (m) Streambed thickness, M (m) Water level in the stream, Hw (m)
10,000.0 50.0 1000.0 0.1 0.05 8 100.0 10.0 0.0 0.5 2.5 0.5 0.5
aquifer connection, based on estimated drawdown and stream depletion. The drawdown and stream depletion curves for HuntÕs solution are not strictly valid after the time when the stream becomes perched above the water table but are shown to illustrate the differences in the two solutions. Fig. 7 shows that drawdown at a location 10 m from the stream is greater when accounting for unsaturated hyporheic zone flow compared to HuntÕs [19] equation. Stream depletion is shown to be less. Unsaturated gravity-driven flow begins at a location on the stream nearest the pumping well at approximately 2.5 days after the initiation of pumping. The maximum specific discharge from the unsaturated segments, jqmax j, is 1.08 m/day and the ultimate interface capillary pressure head, hcu , is 0.08 m. In HuntÕs solution, the specific discharge is 1.0 m/day when the water table is just at the bottom of the streambed. However, assuming flow stays saturated beneath the streambed, specific discharge continues to increase with drawdown. The stream length with unsaturated hyporheic zone flow, yhd , is shown in Fig. 8. Fig. 8 shows that the perched stream length increases rapidly after unsaturated flow begins at 2.5 days. This curve can be explained because the total contributing length, as defined by an effective radius of influence along the stream, increases with time. As the stream becomes perched, stream discharge in the perched length becomes limited to the maximum discharge, qmax . Because the well cannot draw more water from the stream in the nearest reach, it draws more water from aquifer storage, corresponding to larger drawdown compared to the saturated flow case where specific discharge is not limiting. As drawdown in the aquifer increases, the radius of influence increases. Thus, a limited constant specific discharge in the per-
Fig. 7. Comparison between the saturated/unsaturated hyporheic zone flow equations and Hunt [19] for (a) drawdown, sw , at ðx=L; y=LÞ ¼ ð0:1; 0Þ and (b) stream depletion as a function of discharge rate of the pumping well, Qs =Q for hypothetical stream/aquifer system.
Fig. 8. Half-length of stream with unsaturated hyporheic zone flow, yhd , versus time for hypothetical stream/aquifer system.
998
G.A. Fox, D.S. Durnford / Advances in Water Resources 26 (2003) 989–1000
ched reach results in more of the stream contributing to well discharge compared to the saturated case. An intriguing question related to Figs. 7 and 8 is what happens when other wells are pumping from the same aquifer. Additional pumping can create interference with the expansion of the radius of influence. If only a finite stream reach is available to meet the pumped well discharge requirements, the length of stream with an unsaturated hyporheic zone will increase until the entire finite stream is perched. At that time, stream depletion would reach a limit and additional well discharge would need to be provided by release of aquifer storage only, resulting in even larger drawdowns compared to those in Fig. 7. These larger drawdowns have implications for pump and well efficiencies. Cases can exist where the increase in specific discharge between the saturated and unsaturated flow scenarios can result in stream leakage satisfying the flow demand created by pumping. The capillary pressure forces that are created during unsaturated flow will create a greater specific discharge compared to the specific discharge during saturated flow. Therefore, a range of physical conditions exists such that unsaturated flow could satisfy the flow demand created by the groundwater extraction well with less drawdown. The sensitivity of drawdown and stream depletion toto streambed hydraulic conductivity, Ksb , streambed thickness, M, stream width, W , entry pressure head of the aquifer, he , Brooks–Corey parameter, g, and the water level in the stream, Hw , is evaluated using the proposed solution. The parameter values in Table 2 are varied by ±50% of their initial values. Lower streambed conductivity results in greater drawdown and earlier initiation of unsaturated hyporheic zone flow. Higher streambed conductivity results in larger maximum specific discharge, as can be observed by inspection of Eq. (11), supplying more water to satisfy the aquifer stress created by pumping, and delaying the time at which unsaturated flow begins. Stream depletion is also sensitive to variations in the streambed hydraulic conductivity. The maximum specific discharge is inversely related to the streambed thickness. Therefore, as thickness increases, the specific discharge decreases and drawdown in the aquifer increases, indicating that more well discharge is supplied by aquifer storage. Streambed thickness also affects the time at which unsaturated flow begins. Since drawdown is measured with respect to the water level in the stream, a larger drawdown is needed to induce unsaturated flow for a thicker streambed. The stream width also influences the total amount of water available to recharge the aquifer. Total stream depletion increases proportionately to stream width, all other factors being the same, resulting in less aquifer drawdown needed to satisfy the well discharge requirements. This, in turn, delays the onset of unsaturated flow below the stream. The water level in the stream
influences the time when unsaturated flow begins since drawdown is measured with respect to the water level in the stream. A larger drawdown is needed to induce unsaturated flow for deeper streams. The water level in the stream also affects the specific discharge, as shown by Eq. (11). Drawdown and specific discharge were found to be insensitive to the entry pressure head and the Brooks–Corey parameter, g, for the range of parameter values investigated. When streams transition between gaining and losing and/or where episodic pumping leads to frequent transient unsaturated conditions, aeration of the hyporheic zone can occur and influence biogeochemical processes. Dissolved oxygen concentrations will increase, influencing chemical transformations. For example, if pathways exist for air to enter the hyporheic zone below a streambed that has become unsaturated, dissolved oxygen concentrations increase, preventing the development of anaerobic conditions that are beneficial for many nutrient transformations. Anaerobic processes have been shown to be important for organic matter oxidation in some streams [2]. Anaerobic conditions are also important to promote nitrate as an electron acceptor instead of oxygen [13]. Aeration of the hyporheic zone could significantly influence the effectiveness of bank filtration systems for nutrient transformations. Small municipalities across the United States are experiencing increased pressures to improve filtration of surface water resources and further reduce nitrate concentrations in drinking water sources. Instead of constructing improved water quality treatment facilities, many municipalities have investigated the possibility of modifying the treatment process by obtaining source water through groundwater extraction wells located in the alluvium, i.e. a bank filtration process, rather than directly from the surface water supply. In bank filtration, groundwater extraction wells cause water to be interchanged between a stream and alluvial aquifer through a streambed of contrasting hydraulic conductivity. Bank filtration attempts to take advantage of the natural filtration as water flows through the porous material and possible anaerobic conditions that develop in the hyporheic zone beneath the streambed [12]. However, if groundwater extraction leads to unsaturated flow conditions, the aeration of the hyporheic zone could potentially attenuate the aerobic conditions and prevent nutrient transformations.
5. Summary and conclusions The objective of this paper is to discuss the development and implications of unsaturated hyporheic zone flow beneath a stream. Saturated flow is typically assumed for seepage between a stream and an underlying aquifer. However, this flow can become unsaturated if
G.A. Fox, D.S. Durnford / Advances in Water Resources 26 (2003) 989–1000
the water table falls a sufficient distance below a less permeable streambed. Flow from a stream to an underlying aquifer can be characterized by three hydrologic states or regimes: saturated flow, a transition regime, and unsaturated, gravity-driven flow below the streambed. It is shown that the development of a unit hydraulic gradient in the hyporheic zone is rapid and the transition regime can be ignored in analyzing stream/ aquifer interaction. It is also shown that one of the most important differences between saturated and unsaturated flow beneath a streambed is that unsaturated flow transforms streams from constant head boundaries to constant flux boundaries. Equations are developed to calculate the maximum, limiting flux that can occur under unsaturated conditions. The effects of an unsaturated hyporheic zone are illustrated for the case of stream leakage induced by a well pumping from an aquifer that is hydraulically interacting with a partially penetrating stream. It is often assumed that the stream is a constant head boundary and flow between the stream and aquifer is saturated. However, when pumping wells are located sufficiently close to a stream, gravity-driven unsaturated flow can occur and the stream depletion rate is no longer a function of the drawdown. An approximating solution for drawdown and stream depletion is presented for transient, unsaturated hyporheic zone flow. Linear superposition is used to derive equations for aquifer drawdown and stream depletion when portions of the stream are perched. The proposed equations for aquifer drawdown and stream depletion are illustrated for a hypothetical stream/aquifer system and results are compared to an analytical solution that assumes saturated flow. Unsaturated hyporheic zone flow limits the specific discharge of water from the stream to satisfy aquifer stress. Perhaps more importantly, the length of stream that becomes perched can be predicted. It is shown that unsaturated hyporheic zone flow in response to a pumping well results in a larger contributing stream length and, in most cases, more aquifer drawdown when compared to a solution that assumes the stream/aquifer connection remains saturated. Cases can exist where the increase in specific discharge between the saturated and unsaturated flow scenarios can result in stream leakage satisfying the flow demand created by pumping.
Acknowledgements The authors acknowledge the support of the USDA National Needs Fellowship Program, Grant no. 0038420-8825, for supporting this research. The authors also acknowledge the support of the Northern Colorado Water Conservancy District. The authors thank two anonymous reviewers for suggestions that greatly improved the final manuscript.
999
References [1] Allan JD. Stream ecology: structure and function of running waters. Boston, MA: Kluwer Academic Publishers; 1995. [2] Baker MA, Dahm CN, Valett HM. Anoxia, anaerobic metabolism, and biogeochemistry of the stream-water–ground-water interface. In: Jones JB, Mulholland PJ, editors. Streams and ground waters. San Diego, CA: Academic Press; 2000. p. 260–80. [3] Bear J. Dynamics of fluids in porous media. New York: American Elsevier Publishing Company, Inc; 1972. [4] Bouwer H. Groundwater hydrology. New York: McGraw-Hill Book Company; 1978. [5] Brooks RH, Corey AT. Hydraulic properties of porous media. Colorado State University Hydrology Paper no. 3, March 1964. p. 27. [6] Butler JJ, Zlotnik VA, Tsou MS. Drawdown and stream depletion produced by pumping in the vicinity of a partially penetrating stream. Ground Water 2001;39(5):651–9. [7] Calver A. Riverbed permeabilities: information from pooled data. Ground Water 2001;39(4):546–53. [8] Carsel RF, Parrish RS. Developing joint probability distributions of soil water retention characteristics. Water Resour Res 1988; 24(5):755–69. [9] Conrad LP, Beljin MS. Evaluation of an induced infiltration model as applied to glacial aquifer systems. Water Resour Bull 1996;32(6):1209–20. [10] Corey AT. Mechanics of immiscible fluids in porous media. Highlands Ranch, CO: Water Resources Publications; 1994. [11] Corey AT. Pore-size distribution. In: van Genuchten MT, Leij FJ, Lund LJ, editors. Indirect methods for estimating the hydraulic properties of unsaturated soils. California: Riverside; 1992. p. 37– 44. [12] Doussan C, Ledoux E, Detay M. River-groundwater exchanges, bank filtration, and groundwater quality: ammonium behavior. J Environ Qual 1998;27:1418–27. [13] Duff JH, Triska FJ. Nitrogen biogeochemistry and surfacesubsurface exchange in streams. In: Jones JB, Mulholland PJ, editors. Streams and ground waters. San Diego, CA: Academic Press; 2000. p. 197–217. [14] Fox GA. An improved unsaturated flow module for MODFLOWÕs RIVER package. In: Ramirez JA, editor. Proceedings of the 23rd Annual American Geophysical Union Hydrology Days, Colorado State University, Fort Collins, Colorado, 2003. [15] Fox GA, DuChateau P, Durnford DS. Analytical model for aquifer response incorporating distributed stream leakage. Ground Water 2002;40(4):378–84. [16] Glover RE. Mathematical derivations as pertain to groundwater recharge. Fort Collins, CO: Agricultural Research Service, USDA; 1960. [17] Glover RE, Balmer CG. River depletion from pumping a well near a river. Am Geophys Union Trans 1954;35(3):468–70. [18] Hantush MS. Wells near streams with semipervious beds. J Geophys Res 1965;70(12):2829–38. [19] Hunt B. Unsteady stream depletion from ground water pumping. Ground Water 1999;37(1):98–102. [20] Hunt B, Weir J, Clausen B. A stream depletion field experiment. Ground Water 2001;39(2):283–9. [21] Larkin RG, Sharp JM. On the relationship between river-basin geomorphology, aquifer hydraulics, and ground-water flow direction in alluvial aquifers. Geol Soc Am Bull 1992;104:1608–20. [22] Rawls WJ, Brakensiek DJ. Estimation soil water retention from soil properties. J Irrigat Drainage Div, ASCE 1982;108(IR2):166– 77. [23] Rosenshein JS. Hydrology of North America, Region 18, Alluvial valleys. In: Back W, Rosenshein JS, Seaber PR, editors. The Geology of North America, O-2. Boulder, CO: Geological Society of North America; 1988. p. 165–7.
1000
G.A. Fox, D.S. Durnford / Advances in Water Resources 26 (2003) 989–1000
[24] Rushton K. Discussion of ‘‘Unsteady stream depletion from ground water pumping’’ by B. Hunt. Ground Water 1999;37(6): 805. [25] Sophocleous M, Koussis A, Martin JL, Perkins SP. Evaluation of simplified stream-aquifer depletion models for water rights administration. Ground Water 1995;33(4):579–88. [26] Theis CV. The effect of a well on the flow of a nearby stream. Am Geophys Union Trans 1941;22(3):734–8. [27] Theis CV. The relationship between the lowering of the piezometric surface and the rate and duration of discharge from a well
using ground-water storage. Trans Am Geophys Union 1935;2: 519–24. [28] van Genuchten MT. A closed-form equation for predicting the hydraulic conductivity of unsaturated soil. Soil Sci Soc Am J 1980;44:892–8. [29] Winter TC, Harvey JW, Franke OL, Alley WM. Ground water and surface water: a single resource. US Geol Survey Circular 1988;1139:79. [30] White NF, Sunada DK, Duke HR, Corey AT. Boundary effects in desaturation of porous media. Soil Sci 1972;113(1):7–12.