Next Article in Journal
Uncertainty Quantification in Water Resource Systems Modeling: Case Studies from India
Previous Article in Journal
Numerical Analysis of Drag Force Acting on 2D Cylinder Immersed in Accelerated Flow
 
 
Font Type:
Arial Georgia Verdana
Font Size:
Aa Aa Aa
Line Spacing:
Column Width:
Background:
Article

Groundwater–Surface Water Interaction—Analytical Approach

by
Marek Nawalany
1,†,
Grzegorz Sinicyn
1,
Maria Grodzka-Łukaszewska
1,* and
Dorota Mirosław-Świątek
2
1
Faculty of Building Services, Hydro and Environmental Engineering Warsaw University of Technology, 00-653 Warsaw, Poland
2
Institute of Environmental Engineering, Warsaw University of Life Sciences–SGGW, 02-787 Warsaw, Poland
*
Author to whom correspondence should be addressed.
In the memory of Professor Marek Nawalany (1947–2020).
Water 2020, 12(6), 1792; https://doi.org/10.3390/w12061792
Submission received: 7 May 2020 / Revised: 15 June 2020 / Accepted: 19 June 2020 / Published: 23 June 2020
(This article belongs to the Section Hydrology)

Abstract

:
Modelling of water flow in the hyporheic zone and calculations of water exchange between groundwater and surface waters are important issues in modern environmental research. The article presents the Analytical Hyporheic Flux approach (AHF) permitting calculation of the amount of water exchange in the hyporheic zone, including vertical water seepage through the streambed and horizontal seepage through river banks. The outcome of the model, namely water fluxes, is compared with the corresponding results from the numerical model SEEP2D and simple Darcy-type model. The errors of the AHF model, in a range of 11–16%, depend on the aspect ratio of water depth to river width, and the direction of the river–aquifer water exchange, i.e., drainage or infiltration. The AHF model errors are significantly lower compared to the often-used model based on vertical water seepage through the streambed described by Darcy’s law.

1. Introduction

Surface waters and groundwater are elements of the environment that are not isolated from each other. Continuous water exchange with varying intensity occurs between them. The phenomenon can be investigated in the physical context and as an important constituent of water management decision-making. Water exchange within river–aquifer systems can be assessed either numerically or analytically, depending on the feasibility of the applied approach for the assumed purpose. This article presents both approaches with particular emphasis on analytical modelling of the phenomenon. The expected advantage of using analytical models lies in their accuracy and computational simplicity. Stream–subsurface water exchange is currently recognised as a fundamental process affecting the transport and fate of contaminants and other ecologically relevant substances in streams [1]. The quantitative and qualitative characteristics of such exchange are inextricably influenced by several physiographic, climatic, and anthropogenic factors [2,3,4,5,6,7,8,9,10].
A great deal of environmental issues concerning water exchange between rivers and groundwater are resolved by means of various numerical models depending on the assumed spatial scales of flow phenomena. The literature on the subject provides examples of the application of 3D models [11] and analyses based on simulations with two-dimensional models [12]. Despite an increase in the ability of such models to reproduce the complexity of the described processes, water exchange in the hyporheic zone is still subject to insufficiently detailed investigation [13,14]. The river–aquifer interaction is an example where accurate assessment of all three components of groundwater velocity beneath and in the vicinity of a riverbed is essential for correct calculation of flow paths and discrimination of bottom and bank water exchange in the hyporheic zone [14,15]. Because mathematical description of the physical interaction between surface and subsurface flows is a relatively difficult task, a simple approach is adopted in a multitude of practical cases. The two types of flow are modelled separately, and their coupling is accomplished by means of the linear Darcy-type formula accounting for the difference in height of water tables in the river and the adjacent aquifer, and assuming inverse proportionality to the resistance of the riverbed sediment layer separating the two flow domains [14,16]. The formula imitates vertical seepage through the semi-pervious layer, as proposed already in 1979 by Rushton and Tomlinson [17]. The majority of numerical models of regional groundwater flow currently still describe water exchange within river–aquifer systems as vertical water seepage through riverbed sediments, whereas the corresponding water flux (per unit length of river stretch segment) is approximated by the Darcy-type model (DM) represented by Equation (1):
Q t o t = W r ˜ · k s d s ( Φ * H r ) ,   ( m 3 / s / m )
where W r ˜ is the width of the riverbed (m); ks is the hydraulic conductivity of river sediments (m/s); ds is the thickness of river sediments under the river bottom (m); Φ* is the piezometric head in the aquifer next to the riverbed (m); and Hr is the height of water table in the river (m).
This approach is implemented in the widely used MODFLOW model providing a stable and convergent numerical solution. In DM model, Equation (1), friction that hampers water flow through river sediments is usually described by just one lumped parameter—resistance c = d s k s (s). Although the flow intensity through the sediments depends on local minor changes in the values of their hydraulic conductivity, in most groundwater analyses, it is impractical to quantify such small-scale variability, and the seepage Equation (1) is commonly considered a satisfactory basis to approximate the bulk flow rate of water through river sediments [17].
Physically, groundwater flows in three dimensions of space. Models describing water exchange process, however, are frequently reduced to two-dimensional flow in a vertical plane perpendicular to riverbed. The legitimacy of the approach was confirmed by Rushton and Rathod [18]. Technically, any 2D approximation to flow made for a central vertical plane of a homogenous segment of a river stretch is replicated unchanged along each segment. In the majority of regional hydrogeological studies such an approach assumes that the exchange of water between the river and the aquifer occurs only through the river bottom. This simplification is justified when the dimensions of the river banks are significantly smaller when compared to the dimensions of the river bottom. When the size of the river banks is comparable with the width of the riverbed, however, the share of water flow through river banks in total water exchange also increases. Therefore, the simplified estimate of the total water flow through river sediments based on Equation (1) may be significantly in error. Errors in water exchange within the river–aquifer systems are addressed in [19], examining the issue of consistency of water balances of surface waters and groundwater. In this article, analytical and numerical models of water exchange through the river bottom and banks are considered as an essential part of multistage research aimed at modelling of water fluxes in the hyporheic zone of rivers. The research particularly aims at evidencing that the application of analytical models of flow across the river banks in addition to flow through the river bottom may reduce errors in assessing water exchange within river–aquifer systems. This is particularly evident in juxtaposition of the outcomes of the analytical model (or its numerical counterpart) to the measure of water exchange calculated by means of the third type of models—the Darcy-like Equation (1)—when the latter is applied in physically unjustifiable situations.
The analytical approach addressed in this article applies the “method of fragments” developed by Carslaw and Jeager [20]. The method was originally applied to solve analytical problems of steady 2D heat conduction in solid bodies of complex shape [21]. The method was soon undertaken by Polubarinova [22,23]. She adopted and extended the methods of Carslaw and Jeager for solving problems of flow within complex porous structures (aquifers, dikes, etc.). The analytical description of water exchange between the river and aquifer has been reviewed from a historical point of view in the available literature [24]. Similar solutions for the improvement of the accuracy of the river–groundwater system interaction have also been published [25,26].
Through proposing the Analytical Hyporheic Flux model (AHF), this article focuses on an analytical description of the local river–riparian aquifer water exchange, and challenges the insufficient accuracy of the commonly used Equation (1). The analytical formulae of the AHF model permit estimating water fluxes including vertical water seepage through the streambed as well as seepage through river banks.

2. Materials and Methods

2.1. The Flow System

Water exchange between the river and aquifer occurs both through the river bottom and river banks. For rivers where Lbott >> Lbank (Figure 1a), vertical water seepage is the dominant part of total flow. For rivers where both lengths are comparable (Figure 1b), both streams are significant. Flow through the river banks should be particularly included in the description of the river–riparian aquifer water exchange in small rivers of the lowland landscape (e.g., the Upper Biebrza River in Poland) characterised by high depth relative to river width. The latter geometry is typical for small lowland rivers in agricultural regions in Poland. In the developed AHF model, allowing for calculation of the amount of water exchange in the hyporheic zone including vertical water seepage through the streambed and horizontal seepage through river banks, the geometry of the river cross-section and shape of river sediment is approximate with rectangle geometry. The analytical model of water exchange was developed with the assumption of simplified rectangular geometry of the riverbed cross-section (Figure 2).
Figure 3 illustrates the simplified geometry of the river–aquifer system and its relevant variables and parameters. For the sake of simplicity, the flow domain was subdivided into four subdomains, namely P.1, P.2, P.3, and P.4, each with a rectangular shape. Although the partition of the flow domain into simpler fragments is a relatively old method [20,27], it is still used in numerous cases to simplify the analytical approach [28], or in benchmarking of numerical models [29].
In the exemplary system, a vertically averaged (therefore constant) piezometric head in the riparian aquifer Φ * along the right vertical boundaries of subdomains P.1 and P.3, and the height of water table in the river– H r are considered as external enforcing factor for water flow in the river–aquifer system.
The following dimensions are ascribed to river sediments of the river–aquifer system presented in Figure 3: Da is the thickness of the aquifer beneath river sediments (m), Wrs is the half-width of river sediments (m), ds is the thickness of river sediments below the river bottom (m), Wr is the half-width of the riverbed (m), Lar is the one-side length of the riparian aquifer (m), Dar is the maximum thickness of the riparian aquifer (m). The mineral bedrock beneath the riparian aquifer is adopted as the reference level for all variables. In the analytical model of groundwater flow in and under sediments, two state variables are defined: Φa(y) is the piezometric head in the semi-confined aquifer under the river sediments (m), h(y) is the height of free water table in river sediments of P.3 (m), hs is the height of free water table in P.3 at the river bank (m). Aquifer and river sediments are also described in terms of their hydraulic parameters: k a is the hydraulic conductivity of aquifer under river sediments, i.e., in P.1 and P.2 (m/s), k s is the hydraulic conductivity of river sediments in P.3 and P.4 (m/s), c = d s k s is the resistance to water flow in P.4 (s). Moreover, two variables describe external enforcing factors of the river is the aquifer system is the Figure 3: Φ* is the averaged piezometric head in the riparian aquifer along the right vertical boundary of subdomains P.1 and P.3 (m), H r is the height of water table in the river (m). Two resultant variables are calculated by the AHF model: Q b o t t is the vertical seepage through the river bottom per 1 meter of the river stretch segment from/to flow domain P.2 (m3/s/m), Q b a n k is the horizontal seepage through the river bank per 1 meter of the river stretch segment from/to flow domain P.3 (m3/s/m).
Two halves of the river sedimentary envelope and two halves of a riparian valley—left (L) and right (R)—are neither physically nor geometrically identical. The model of seepage in river sediments and groundwater flow in the adjacent aquifer generally requires assessment of two different sets of parameters corresponding to L and R parts of the river–aquifer system. The assessment of the shape and location of a no-flow boundary between L and R parts of the system, however, is not possible with a simple analytical model considered in this article, because it requires solving a free boundary problem instead. In this article, river–aquifer water exchange is assessed by formulating an analytical model with only one set of identical parameters for the L and R sides of the flow system. Future stages of the research are planned to extend the present analytical model through relaxation of its most restrictive assumptions. The model is therefore anticipated to become suitable for more accurate assessment of water exchange with better approximation to the real flow system.

2.2. The AHF Model–Analytical Description of Water Seepage through River Sediments

Groundwater flow equations and their solutions used for calculating water exchange within different types of river–aquifer systems have been presented in the literature for a relatively long time. Water exchange through the river bottom has been in the focus of analytical models in a multitude of cases. A distinctive feature of the AHF model is simultaneous consideration of water exchange through sediments below the river bottom–P.4, and through sediments to the right of the river bank–P.3 (Figure 3). Several assumptions were made in order to simplify the analytical model of seepage in river sediments in terms of its physical adequacy and computational feasibility. The analytical model is aimed to ensure that the continuity of flow between the two water environments is fulfilled. The assumptions refer to geometry and processes in subdomains P.1 and P.2 of the aquifer under river sediments, and P.3, P.4 in river sediments. The simplifying assumptions for the AHF analytical model are as follows:
  • Due to its slow dynamics, groundwater flow in and under river sediments can be approximated with a steady state flow. The assumption was analysed by Nawalany [30] who evidenced that despite a rapidly fluctuating river water table, a slow response of groundwater flow in the adjacent aquifer occurs as a consequence of dumping of high frequency changes in pore pressure by porous rocks. The adequacy of such an assumption was also discussed in more recent literature [26,31,32]. Preliminary field measurements in a real river–aquifer system show that water table fluctuations in the river– H r , and slow response of free water table that follows in the adjacent riparian aquifer– Φ * support the assumption of merely quasi-steady water exchange in the exemplary river–aquifer system. Both variables H r and Φ * are used in the analytical AHF model as two independent external variables enforcing water flow in river sediments.
  • The symmetry of the L and R parts of the river–aquifer system implies an existence of a no-flow boundary within the aquifer bellow river sediments. In the AHF model, this boundary is assumed to be a vertical line in the middle of the river. This rather restrictive simplification will be relaxed in the future developments of the model.
  • Flow in subregion P.4 is assumed to be approximately vertical, because the left and the right boundaries of P.4 are vertical, whereas the upper boundary condition of the river water table is horizontal. Subregion P.4 is therefore described as a semi-pervious layer exerting a resistance drag (c) to water flowing from P.2 to the river bottom.
  • Water flow between P.3 and P.4 is assessed as negligible. Therefore, the internal boundary between the two subregions is considered a no-flow boundary. The boundary between P.1 and P.3 is also assumed to be a no-flow boundary.
  • The base of the river–aquifer system is impervious. Direct recharge of river sediments from the top by infiltrating precipitation is also assumed negligible.
Below flow in P.1, P.2, and P.3 parts of the river–aquifer system is described in terms of their boundary conditions, groundwater flow equations and their solutions, and respective parameters. Flow in subdomain P.4 is approximated by vertical seepage through the semi-pervious layer, and is linked to flow in P.2.
Flow in subdomain P.1 is confined and described by the Laplace equation:
k a ( 2 Φ a ( y , z ) y 2 + 2 Φ a ( y , z ) z 2 ) = 0 ,   for   y [ W r ,   W r s ] ,   z [ 0 ,   D a ]
with boundary conditions: left boundary is the horizontal component of specific discharge k a   Φ a ( W r + ,   z   ) y equal to the horizontal component of specific discharge k a   Φ a ( W r ,   z   ) y at the vertical boundary between P.1 and P.2 and piezometric head continuous, i.e.,
Φ a ( W r ,   z   ) = Φ a ( W r + ,   z )
right boundary is the imposed piezometric head (1st type b.c.), i.e.,
Φ a ( W r s , z ) = Φ * = c o n s t
lower boundary is the impervious boundary (2nd type b.c.), i.e.,
Φ a ( y , 0 ) z = 0
upper boundary is the no flow boundary (2nd type b.c.), i.e.,
Φ a ( y ,   D a ) z = 0
because no water exchange is assumed between P.3 and P.1.
The solution to Equation (2) can be found by means of the factorisation method like in Nawalany‘s work [27], where some non-zero penetration of the riverbed into a confined aquifer D r > 0 was assumed. When D r converges to zero (as in the example of this article), the solution converges to
Φ ˜ a ( y , z ) = i = 1 D i sinh [ μ i ( W r s y ) ] cos ( μ i z ) + Φ ˜ * β * ( W r s y ) for   y [ W r ,   W r s ] ,   z [ 0 ,   D a ]
where
Φ ˜ a ( y , z ) = Φ a ( y , z ) H r
Φ ˜ * = Φ * H r
whereas its derivatives are equal to
Φ ˜ a ( y , z ) y = i = 1 D i μ i cosh [ μ i ( W r s y ) ] cos ( μ i z ) + β *
Φ ˜ a ( y , z ) z = k = 1 D i μ i sinh [ μ i ( W r s y ) ] sin ( μ i z ) .
Elements of series ( D i ) ,   i = 1 , 2 . and value of β * are derived further from requirements of flow continuity at the boundary between P.1 and P.2.
Solution of Equation (7) satisfies right b.c. as Φ ˜ a ( y = W r s , z ) = Φ ˜ * and, consistently, Φ ˜ a ( y = W r s , z ) z = 0 for all z [ 0 ,   D a ] . For y [ W r ,   W r s ] , Equation (9) satisfies lower b.c., i.e., Φ ˜ a ( y ,   z = 0 ) z = 0 . To also satisfy the no flow condition on the upper boundary of P.1, i.e., Φ ˜ a ( y ,   z = D a ) z = 0 , μ i must be equal to
μ i = i   π D a   for   i = 1 , 2 , ,
Notice also that total outflow from P.1 into the riparian aquifer is equal to
Q P . 1 r i p = k a 0 D a Φ ˜ a ( y = W r s , z ) y d z = k a i = 1 D i μ i 0 D a cos ( μ i z ) d z k a D a β * = k a D a β *
The minus sign in Equation (11) means that at that border, water physically flows along the x-axis if β * < 0 , and opposite to y-axis if β * > 0 .
Flow in subdomain P.2, i.e., in the part of the aquifer under river bottom sediments of P.4, is also described by the Laplace equation:
k a ( 2 Φ a ( y , z ) y 2 + 2 Φ a ( y , z ) z 2 ) = 0   for   y [ 0 ,   W r ] ,   z [ 0 ,   D a ]
with boundary conditions: left boundary is the no flow boundary (2nd type b.c.), i.e.,
Φ a ( 0 , z ) y = 0
right boundary is the horizontal component of specific discharge k a   Φ a ( W r ,   z   ) y equal to horizontal component of specific discharge k a Φ a ( W r + ,   z ) y at the vertical boundary between P.2 and P.1, and also piezometric head is continuous, i.e.,
Φ a ( W r ,   z   ) = Φ a ( W r + ,   z )
lower boundary is the impervious boundary (2nd type b.c.), i.e.,
Φ a ( y , 0 ) z = 0
upper boundary is the seepage through a layer of river sediments of P.4 to river bottom, i.e.,
q s = k a Φ a ( y ,   D a ) z = Φ a ( y ,   D a ) H r c
The solution to the flow Equation (12) in P.2 satisfying b.c.-s, Equations (13)–(16), can be found by means of the factorisation method. It is given by the following formula
Φ a ( y , z ) = H r + k = 1 A k cosh ( λ k y ) cos ( λ k z ) ,   for   y [ 0 ,   W r ] ,   z [ 0 ,   D a ]
whereas its derivatives are equal to
Φ a ( y ,   z ) y = k = 1 A k λ k sinh ( λ k y ) cos ( λ k z )
Φ a ( y ,   z ) z = k = 1 A k λ k cosh ( λ k y ) sin ( λ k z )
Parameters λ k are specified through b.c. (16)
k a Φ a ( y ,   D a ) z = Φ a ( y ,   D a ) H r c   Φ ˜ a ( y , D a ) z = Φ ˜ a ( y , D a ) χ 2
where   Φ ˜ a ( y , z ) = Φ a ( y , z ) H r and   χ 2 = k a   c
Substituting Equation (17) and Equation (19) to Equation (20) provides k = 1 A k λ k χ 2 cosh ( λ k y ) tg ( λ k D a ) cos ( λ k D a ) = k = 1 A k cosh ( λ k y ) cos ( λ k D a ) = > λ k χ 2 tg   ( λ k D a ) = 1 Elements of series ( λ k ) ,   k = 1 , can be calculated from the non-linear algebraic equation
D a χ 2 1 λ ˜ k = t g ( λ ˜ k )
where   λ ˜ k = λ k D a
for instance by the Newton method.
Elements of series ( A k ) ,   k = 1 , 2 , are derived below from the requirements of continuity along the vertical boundary between P.1 and P.2. Once parameters A k and λ k are known,   Q b o t t can be calculated from Equation (19) by integrating specific discharge q z over top horizontal boundary of subregion P.2, i.e.,
Q b o t t = 0 W r q z ( y , D a ) d y = 0 W r k a Φ a ( y ,   D a ) z d y = k a k = 1 A k λ k sin ( λ k D a ) 0 W r cosh ( λ k y ) d y
hence
Q b o t t = k a k = 1 A k sin ( λ k D a ) sinh ( λ k W r ) ,   ( m 3 / s / m )
Flow continuity between subregion P.1 and P.2. In order to evaluate elements of series ( A k ) ,   k = 1 , 2 , , ( D i ) ,   i = 1 , 2 , and value of β * , the piezometric head and horizontal components of specific discharge must be assumed continuous along the boundary between flow subregions P.1 and P.2, i.e., at   y = W r and for z [ 0 ,   D a ]
k = 1 A k cosh ( λ k W r ) cos ( λ k z ) = i = 1 D i sinh [ μ i ( W r s W r ) ] cos ( μ i z ) + Φ ˜ * β * ( W r s W r )
k = 1 A k λ k sinh ( λ k W r ) cos ( λ k z ) = i = 1 D i μ i cosh [ μ i ( W r s W r ) ] cos ( μ i z ) + β *
Once elements of series ( A k ) ,   k = 1 , 2 …are known, value of β * can be calculated from mass conservation in flow subregions P.1 and P.2, i.e., from equality of Equations (11) and (25)
Q P . 1 rip = Q bott k a D a β * = k a k = 1 A k sin ( λ k D a ) sinh ( λ k W r ) β * = k = 1 A k sin ( λ k D a ) sinh ( λ k W r ) / D a
To evaluate elements of the two series, ( A k ), k = 1 , 2 … and ( D i ) , i = 1 , 2 , equality (6.2) needs to be multiplied by b = ( W r s W r ) and added to Equation (26), resulting in
k = 1 A k [ cosh ( λ k W r ) + λ k b sinh ( λ k W r ) ] cos ( λ k z ) = = j = 1 j = D j [ sinh ( μ j b ) μ j b   cos h ( μ j b ) ] cos ( μ j z ) + Φ ˜ *
Ultimately, values of ( 2 N + 2 ) unknowns— ( A k ), k = 1 , 2 , ( N + 1 ) , ( D i ) , i = 1 , 2 , ( N + 1 ) —need to be found. Simultaneous substitution of Equation (28) for β * and z m = m D a N ,   m = 0 , 1 ,   N for z into Equations (26) and (27) leads to a set of ( 2 N + 2 ) linear algebraic equations, where infinite series are approximated with their finite counterparts
k = 1 N + 1 A k { cosh ( λ k W r ) cos ( λ k z m ) b D a   s i n ( λ k D a ) sinh ( λ k W r ) } i = 1 N + 1 D i sinh [ μ i b ] cos ( μ i z m ) = Φ ˜ * , m = 0 , 1 , ,   N
k = 1 N + 1 A k { λ k sinh ( λ k W r ) cos ( λ k z m ) + s i n ( λ k D a ) s i n h ( λ k W r ) / D a } + + i = 1 N + 1 D i μ i cosh [ μ i   b ] cos ( μ i z m ) = 0   , m = 0 , 1 , ,   N
Hyperbolic sine and cosine in Equations (30) and (31) need some pre-calculation to avoid rising exponents, i.e., s i n h   α = e x p   α [ 1 e x p   ( 2 α ) ] / 2 and c o s h   α = exp α [ 1 + e x p   ( 2 α ) ] / 2 . By denoting
A k ^ = A k e x p   ( λ k W r )
D ^ i = D i e x p   ( μ i b )
  • ω ˜ m k = { cos ( λ k z m ) b D a   s i n ( λ k D a ) t g h ( λ k W r ) } [   1 + e x p   ( 2 λ k W r )   ] / 2
  • ω ˜ ˜ m k = { { λ k cos ( λ k z m ) + s i n ( λ k D a ) D a } [   1 e x p ( 2 λ k W r ) ] / 2
  • δ ˜ m i = cos ( μ i z m ) [ 1 e x p ( 2 μ i b ) ] / 2
  • δ ˜ ˜ m i = μ i cos ( μ i z m ) [ 1 + e x p ( 2 μ i b ) ] / 2
Equations (30) and (31) can be written shortly as
k = 1 N + 1 A k ^ ω ˜ m k i = 1 N + 1 D ^ i δ ˜ m i = Φ ˜ * , m = 0 , 1 , ,   N
k = 1 N + 1 A k ^ ω ˜ ˜ m k + i = 1 N + ! D ^ i δ ˜ ˜ m i = 0   , m = 0 , 1 , ,   N
or in the block-matrix form for as
Ω ˜ N + 1 x N + 1 | Ω ˜ ˜ N + 1 x N + 1 |   Δ ˜ N + 1 x N + 1 Δ ˜ ˜ N + 1 x N + 1   A ^ N + 1 x 1 D ^ N + 1 x 1 = Ρ 2 N + 2 x 1
where
  • Ω ˜ ( N + 1 ) x   ( N + 1 ) = { ω ˜ m k } ( m = 0 , 1 , ,   N   ;   k = 1 , , N + 1 )
  • Ω ˜ ˜ ( N + 1 ) x   ( N + 1 ) = { ω ˜ ˜ m k } ( m = 0 , 1 , ,   N   ; k = 1 , , N + 1 )
  • Δ ˜ ( N + 1 ) x   ( N + 1 ) = { δ ˜ m k }   ( m = 0 , 1 , ,   N   ; k = 1 , , N + 1 )
  • Δ ˜ ˜ ( N + 1 ) x   ( N + 1 ) = { δ ˜ ˜ m k } ( m = 0 , 1 , ,   N   ; k = 1 , , N + 1 )
  • A ^ ( N + 1 )   x   1 = [ A k ^ ]   ( k = 1 , , N + 1 ) is the first (sub)vector of unknowns
  • D ^ ( N + 1 )   x   1 = [ D ^ i ] ( i = 1 , , N + 1 ) is the second (sub)vector of unknowns
  • Ρ ( 2 N + 2 )   x 1 = [ Φ ˜ * _   0   _ ] is the right hand side vector.
After solving Equation (36) for ( A k ) ,   k = 1 , 2 , ( N + 1 ) , the sought seepage through river bottom is the Q b o t t can be readily calculated from approximations of Equation (11) or (24)
Q b o t t = k a D a β *   or
Q b o t t = k a k = 1 N + 1 A k sin ( λ k D a ) sinh ( λ k W r ) = k a k = 1 N + 1 A k ^ sin ( λ k D a ) [ 1 e x p   ( 2 λ k W r ) ] / 2
Convergence of finite series ( A k ) ,   k = 1 , 2 , ( N + 1 ) and ( D i ) ,   i = 1 , 2 , ( N + 1 ) has been checked for increasing N .
Flow through river sediments in subdomain P.3 is unconfined and, after moving the origin of coordinates to point ( y o = W r ,   z o = D a ) , described by the Laplace equation
k s ( 2 Φ s ( y , z ) y 2 + 2 Φ s ( y , z ) z 2 ) = 0   ,   for   y [ 0 ,   b ] ,   z [ 0 , h ( y ) ]
where
Φ s ( y , z ) = Φ ( y , z ) D a is the piezometric head in sediments of P.3 (m),
b = W r s W r   is the width of P.3 (m),
h ( y ) : = h ( y ) D a is the height of free water table in P.3 over the reference level z 0 = D a , (m)
with boundary conditions at:
left boundary consisting of three segments: (assumed) no-flow boundary between P.3 and P.4 (2nd type b.c.), i.e., Φ s ( 0 ,   z ) y = 0 ,   0 z d s , constant head along the water part of the river bank (1st type b.c.), i.e., ϕ s ( 0 , z ) = H r , d s < z H r , seepage face h s along the aerial part of the river bank (1st type b.c.), i.e.,
ϕ s ( 0 ,   z ) = z ,   H r < z h s = h ( 0 )
right boundary is the constant piezometric head (1st type b.c.), i.e.,
ϕ s ( b ,   z ) = ϕ * ,   0 z ϕ *
lower boundary is the (assumed) no flow boundary between P.1 and P.3 (2nd type b.c.), i.e.,
  ϕ s ( y ,   0   ) z = 0 ,   0 y b
upper boundary is the free boundary of P.3, i.e., ϕ s ( y , h ( y ) ) = h ( y ) with no recharge from the top i.e.,
  ϕ s ( y ,   h ( y ) ) z = 0 ,   0 y b
External constraints are also referenced to z 0 = D a , i.e., H r : = H r D a and ϕ * : = ϕ * D a .
Because the height of free water table h ( y ) in P.3 is unknown, and so is the height of seepage face h s = h ( 0 ) , algebraic relationship ϕ s ( y , h ( y ) ) = h ( y ) is added to make the solution to Equation (39) unique. The unknown shape of subregion P.3 makes the problem Equations (39)–(43) into a free boundary issue, and is cumbersome to solve. The calculation of river inflow/outflow to/from subregion P.3,   Q b a n k , however, is possible by means of the method originally derived by Czarny [23]. It allows for calculating unconfined flow in sediments of subdomain P.3 without actual solving of the free boundary problem Equations (39)–(43). By applying Leibnitz theorem to Darcy’s Law, the following general formula for total flow in P.3, Q s ( y ) , can be derived
Q s ( y ) = 0 h ( y ) [ k s Φ s ( y , z ) y ] d z = k s y { 0 h ( y ) Φ s ( y , z )   d z [ h ( y ) ] 2 2 } = = k s y {   h ( y ) Φ ˜ ( y ) [ h ( y ) ] 2 2 }
where Φ ˜ ( y ) = 1 h ( y ) 0 h ( y ) Φ s ( y , z )   d z is the average piezometric head in river sediments of P.3 along any arbitrary vertical at y over the assumed reference level at z 0 = D a .
Equation (44) has been originally derived by Czarny [23] to justify the use of the Dupuit parabola when calculating inflow of groundwater to an excavation. In this article, the Czarny method is slightly generalised by considering a no-flow boundary at lower part of the left vertical boundary of P.3 (see Equation (40)).
Because there is neither seepage from P.1 to P.3 nor recharge from the top, the continuity of flow in P.3 requires that Q s ( y ) y   = 0 , and therefore
2 y 2 [ h ( y ) Φ ˜ ( y ) [ h ( y ) ] 2 2 ] = 0
The general solution to this equation can be proposed in a simple form
h ( y ) Φ ˜ ( y ) [ h ( y ) ] 2 2 = A y 2 + B y + C
Substitution of parabola, Equation (46) to Equation (45) results in A = 0 , and therefore
h ( y ) Φ ˜ ( y ) [ h ( y ) ] 2 2 = B y + C
whereas the sought total flow at any y can be calculated from Equation (44), i.e. ,   Q s ( y ) = k s B .
Constants B and C can be assessed from the assumed boundary conditions at the P.3 left and right boundaries, namely for y = 0: h ( 0 ) Φ ˜ ( 0 ) [ h ( 0 ) ] 2 2 = C , where h ( 0 ) = h s , whereas term h ( 0 ) Φ ˜ ( 0 ) at the river bank is equal to
h ( 0 ) Φ ˜ ( 0 ) = 0 h ( 0 ) Φ s ( 0 , z ) d z = 0 d s Φ s ( 0 , z ) d z + d s H r Φ s ( 0 , z ) d z + H r h s Φ s ( 0 , z ) d z .  
The first integral can be derived from the observation that for ϕ * > H r ,piezometric head along the flowline located at lower no-flow horizontal part of P.3 boundary, i.e., at z = 0   ( z 0 = D a ) , and continuing further along the no-flow lower part of the left vertical boundary of P.3 at y = 0 , changes linearly from ϕ * down to H r . Therefore, for 0 z d s   ,   ϕ s ( 0 , z ) = H r + ϕ * H r b + d s ( d s z , ) and hence 0 d s ϕ s ( 0 , z )   d z = H r d s + ϕ * H r b + d s d s 2 2 . Using b.c. Equation (40) for other parts of the left boundary of P.3, the remaining two integrals can be calculated, d s H r ϕ s ( 0 , z )   d z = H r ( H r d s ) and H r h s ϕ s ( 0 , z ) d z = h s 2 H r 2 2 . Adding the three integrals of Equation (48) provides C = h ( 0 ) Φ ˜ ( 0 ) [ h ( 0 ) ] 2 2 = H r 2 2 + Φ * H r b + d s d s 2 2 . For y =   b : h ( b ) ϕ ˜ ( b ) [ h ( b ) ] 2 2 = B b + C ϕ * ϕ * [ ϕ * ] 2 2 = B b + C B = { [ ϕ * ] 2 2 H r 2 2 ϕ * H r b + d s d s 2 2 } b , and therefore one-side
Q b a n k = Q s ( y = 0 ) = k s B = k s 2 b   { [ Φ * ] 2 H r 2 Φ * H r b + d s d s 2 }
For Φ * > H r ,    Q b a n k < 0 , which indicates inflow to the river from P.3.
For H r > Φ * , flow reverses, and there is no seepage face either at the left nor at the right boundary of P.3, i.e.,   h ( 0 ) = H r and h ( b ) = Φ ˜ ( b ) = Φ * . Then, total flow, Equation (44) can be evaluated from two integrals of Equation (48) - 0 d s Φ s ( 0 , z )   d z and d s H r Φ s ( 0 , z )   d z , as
Q b a n k = k s 2 b   { H r 2 [ Φ * ] 2 H r Φ * b + d s d s 2 }
After moving back to the original global of coordinates ( y o = 0 ,   z o = 0 ) , Equations (49) and (50) can be shortly written as:
Q b a n k = k s 2 b { ( H r D a ) 2 ( ϕ * D a ) 2 H r ϕ * b + d s d s 2 }
When d s = 0 , Equation (51) reduces to the Dupuit formula Q b a n k = k s 2 b { ( H r D a ) 2 [ Φ * ] 2 ( ϕ * D a ) 2 } .
Ultimately, through adding formulae of outflow/inflow from/to riparian aquifer to/from a river through its bottom— Equation (37) or (38)—and bank— Equation (51)—the total one side river recharge per 1 m of the river segment can be calculated as follows:
Q t o t a l = Q b o t t o m + Q b a n k ,   ( m 3 / s / m )

2.3. The SEEP2D Model–Numerical Approximation of Water Seepage in the River–Aquifer System

River–aquifer water exchange calculated with the AHF analytical model can be compared with the corresponding water flows approximated by SEEP2D software suitable for modelling a variety of problems involving seepage. The SEEP2D is a 2D finite element, steady state, flow model successfully applied in cross-section (profile) models representing a vertical slice through a flow system such as earth dams or levees. The SEEP2D model is based on the Equation (53) [33].
y ( K y y h y + K y z h z ) + z ( K z z h z + K z y h y ) = 0
where h is the total head (elevation head plus pressure head), K is the hydraulic conductivity tensor.
SEEP2D permits modelling for the following conditions: isotropic and anisotropic soil properties; confined and unconfined flow for cross-section models; saturated/unsaturated flow for unconfined cross-section models; flow simulation in saturated and unsaturated zones; heterogeneous soil conditions. In the model, the Laplace Equation (53) is solved by means of the finite element method (FEM), and the flow domain is represented by a finite element mesh consisting of triangular and quadrilateral elements. In unconfined problems, where the position of the free surface of water is unknown, SEEP2D allows for modelling of either the deformation of the mesh to the phreatic surface, or simulation of flow in both saturated and unsaturated zones. In the first approach, flow occurs only in the saturated zone, and is iteratively finding the location of the phreatic surface. The mesh is deformed or truncated for the upper boundary of the mesh to match the phreatic surface. When unsaturated flow is simulated, hydraulic conductivity is modified using either the linear frontal method or the Van Genuchten method [34]. The following boundary conditions can be used in the model at the node in the mesh: constant head (Dirichlet boundary condition), head equals the elevation (exit face), flow rate. Known flux rate is used as a boundary condition along a sequence of element edges on the perimeter of the flow domain. Exit face boundary conditions are used when simulating unconfined flow problems, and should be added along the face where the free surface is likely to exit the flow domain. The SEEP2D program calculates the head, flow, discharge (Darcian) velocity, and pore pressure at every node in the mesh. In the example of water seepage considered in this article, the flow domain-2D cross-section-a vertical slice through a flow system, (Figure 3) and flow condition–unconfined flow in the saturated zone, permit applying the SEEP2D model with the deforming mesh option.

3. Simulation Results and Discussion

Six scenarios of the difference in the mutual position of the water table in the aquifer and river were selected for analysis. The initial scenario, considered the most common, is scenario 3 and 4 (0.5 m difference between Φ* and Hr for the gaining/losing type of river). Such differences are observed for small and deep lowland rivers (e.g., the Upper Biebrza River in Poland). During the drainage period, the difference between Φ* and Hr is approximately 0.5 m during infiltration periods, the difference reaches maximum 2 m. Other scenarios were selected to test models in a wide range of differences between Hr and Φ* both when the river is of draining and infiltrating type.
Parameter values presented in Table 1 were used for calculating water exchange flows by means of the AHF and SEEP2D models for the exemplary river–aquifer system. These parameters were estimated based on of the real river–aquifer system of the Upper Biebrza River in Poland.
Performing the simulation with the AHF model requires acquiring the variables listed in Table 1. In addition to the knowledge of the hydraulic parameters, the application of SEEP2D software for calculations requires the development of a 2D finite element mesh for the flow domain. In the calculation example, a finite element mesh representing the river–aquifer system flow domain (Figure 3) was generated with mesh generation tools provided in GMS Groundwater Modelling System. In the model, the “no-flow” boundary was applied to the bottom of the aquifer system. Constant head boundary conditions are in locations where the head and river water table are known. The fixed value of groundwater table in the riparian aquifer Φ* (Figure 3) along the left and right perimeter of the flow domain (where nodes elevation is not greater then Φ*), and fixed value of river table along the wetted perimeter in the banks and bottom of the river. When river is of the gaining type (Hr < 27 m) exit face boundary conditions are placed along the river perimeter (left and right river bank above the water table). When the river is of the losing type (Hr > Φ*), exit face boundary conditions are placed along the left and right perimeter of the flow domain in nodes with an elevation greater than Φ*. In the remaining part of the perimeter of the flow domain, the “no-flow” boundary conditions are assumed.
In calculations where SEEP2D reflected the simplifying assumptions made in the AHF model, the boundary conditions are implemented by placing constant head boundary conditions in locations where the head and river water are known, and “no-flow” boundary conditions on the remaining part of the perimeter of the flow domain (the same boundary conditions were applied as in the AHF model). The assumption of no flow between P.3 and P.4 as well as between P.1 and P.3 (Figure 3) was implemented in SEEP2D software by creating a discontinuity in the mesh—an impermeable flow barrier with a width of 0.05 m—and placing “no-flow” boundary conditions in nodes along the outer perimeter of this barrier.
In the case of exit face boundary conditions, if the head at a node on the boundary becomes greater than the node elevation during the iteration process, the head at the node is fixed at the nodal elevation and the node acts as a specified head boundary. In the calculation, the SEEP2D software was used with the deforming mesh option (unconfined flow in the saturated zone). In this case, the boundary conditions are implemented by placing constant head boundary conditions in locations where the head and river water table are known, and iteratively is finding the location of the phreatic surface, and the mesh is deformed or truncated so that the upper boundary of the mesh matches the phreatic surface, placing exit face boundary conditions along the boundary where the phreatic surface is assumed to exit.
The basic finite element mesh in the SEEP2D model consist of 1185 nodes and 2127 triangle elements. It was refined in an area of sediment close to river banks and bottom, with high flow or high gradient in head. The mesh element size along the river cross-section is Δx = 0.25 m. To ensure a suitable mesh size, test calculations were performed for a twice denser mesh (4564 nodes and 8508 elements, Δx = 0.125 m). A test with decreasing density of the mesh (522 nodes, 944 elements, Δx = 0.5 m) was also performed. Calculation was done for all scenarios presented in Table 2. The estimated differences in calculated flow across the river bottom and banks for the base, refine, and decrease finite element mesh are in a range from 0.27% to 0.57%. These test results show that the base finite element mesh in SEEP2D model has a suitable mesh size.
The results of the analytical AHF model in terms of Qbank, Qbott, and Qtot are listed in Table 3 with the corresponding numerical results of SEEP2D (with assumption of no flow between P.3 and P.4 as well as between P.1 and P.3).
Flow across the river bank (Qbank) estimated with the AHF model varies from 28% to 11% with respect to the numerical solution. The error can be explained by the application of an approximate method for assessing flow in subregion P.3. The flow across the river bottom calculated with the AHF model (Qbott) in all scenarios does not differ from the numerical solution by more than 2%. Such good accuracy can be explained by presenting the exact analytical solution, Equation (25) to flow equations in the form of an infinite series expressed in terms of two infinite series—(Ak), k = 1,2,… and (λk), k = 1,2,… Qbott—has been approximated with a finite series, Equation (38) of length N+1, ( Q b o t t ( N ) ) , in which (Ak), k = 1,2,…, N+1 were calculated from a system of linear algebraic equations, Equation (36). Elements of series ( λ k ) ,   k = 1 , 2 , ,   N + 1   solutions to algebraic Equation (22) were calculated by means of the Newton method using parameters presented in Table 1. Values of the first elements in series ( λ k ) ,   k = 1 , 2 , …= are as follows (0.02777; 0.162368; 0.3168768; 0.473060; 0.6296871; 0.7864939; 0.9433914; 1.1003407; 1.2573226; 1.4143261; 1.5713448,….). The convergence of the finite series ( Q b o t t ( N ) ) was checked for increasing N . Figure 4 shows how N affects Q b o t t ( N ) . It was proven through repetitive calculation of formula (38) that Q b o t t ( N ) , converges to the SEEP2D solution already for N = 200, and later, from N = 200 to N = 3000, this convergence holds. The differences between the models for the number N greater than 200, for some scenarios, can be explained by the fact that the numerical solution obtained using the SEEP2D model, although very accurate, is an approximate solution. Qbott differences between the two models depend on the number N and range from 0.04 × 10−6 m3/s/m (for N = 400, scenarios 3 and 4) to 0.2 × 10−6 m3/s/m (for N = 3000, scenarios 1 and 6). The discrepancy range is from 0.14% to 2% respectively, and does not significantly affect further considerations.
The AHF model also allows for assessing specific discharge and hydraulic head in flow subregions P.1, P.2, P.3, and P.4. (Figure 5).
Components of specific discharge were calculated by means of Darcy Law and formulae for y- and z-derivatives of piezometric head: Equations (8) and (9) in P.1, Equations (18) and (19) in P.2, and Equation (19) in P.4 (assuming flow continuity across top horizontal boundary of P.2). Notice (Figure 5) that the AHF model accurately reflects the specific discharge field in subdomains P.1, P.2, and P.4. Due to the generalisation of Equations (48) and (49) of the method proposed by Czarny [23] for approximating total flow, and Dupuit formula for height of free water table, only an approximate value of the horizontal component of specific discharge can be determined in subdomain P.3. This, however, does not significantly affect the total flow in this subdomain.
The AHF model described in this article was also compared with the numerical model SEEP2D describing a real flow system, i.e., where simplifying assumptions (of no flow between P.3–P.4 and P.1–P.3) are rejected. The results obtained using a simple model based on the Darcy formula (DM model) were also analysed. Table 4 and Figure 6 present the results of the analytical AHF model, Q t o t , in comparison with the corresponding results from the SEEP2D model, calculated for the entire domain (for both sides of the river), and Darcy’s law (DM)—Equation (1).
In further analyses, simulation results obtained with the application of the SEEP2D numerical model without simplifying assumptions were assumed as the reference solution.
The results obtained from this SEEP2D model indicate that the assumption in the analytical model of no flow between P.3–P.4 and P.1–P.3 generates errors from 24% for scenario 1 to 26% for scenario 6. The assumption of no flow between P.3–P.4 and P.1–P.3 reduces the actual flow through the river bottom and river bank. In all calculation scenarios with of the SEEP2D model, water does not leave the model along the exit face on river banks.
For all scenarios, the analytical model and estimation based on Equation (1) underestimate the values of water flow in the hyporheic zone (Figure 7). AHF model errors range from 11% to 16%. Error values in river seepage estimated by means of Darcy’s law are much larger and reach values in a range of 40–48%.
In comparison to the DM model, the AHF model describes seepage through river sediments in more detail, because as it not only considers external variables–Hr, Φ*, but also permits calculation of both components of the specific discharge (qy and qz) in subdomain P.2 located under the river bottom.
It is worth emphasising that, in contrast to the SEEP2D model, bottom recharge calculated by means of the AHF and DM models does not change as to its in absolute value when the river changes its type from draining into the infiltrating one (Figure 8).
Differences in the SEEP2D numerical solutions for a drainage and infiltrating river seen in Figure 9 can be explained by changing the position of free water table, and hence the flow domain shape in both cases. In this situation, the absolute Qbott value decreases by 7% from 29.5 × 10−6 m3/s/m for a drainage river to 27.5 × 10−6 m3/s/m for an infiltrating river.
The AHF model inaccuracy is associated with the adopted simplifications and assumptions. Although the “method of fragments” enabled derivation of consistent analytical solutions to groundwater flow equations for the river–aquifer system, assumptions of no water exchange between flow subdomains P.1 and P.3, and between P.3 and P.4 (Figure 3) have an unavoidable impact on the analytical description of the resultant flow.
Considering seepage through the river bank in the calculations permits a satisfactory approximation of water flow through river sediments in the case of small riverbeds close to a rectangular shape and significant water depths.
In the case of rivers where water depth is comparable to river width, the share of flow through the river bank is a considerable part of total water exchange between the river and sediment layer. In the analysed case, for a depth of 3.5 m, in both AHF and SEEP2D models, Qbank accounts for approximately 40% of total flow (Figure 10).
Therefore, the application of a formula-based model resulting from Darcy’s law (DM model) leads to major errors in this case. In particular for a water depth of 3.5 m, the error of the DM model is 48%.
This type of approach (DM model), implemented e.g., in the MODFLOW, MODBRANCH, MIKE-SHE, and HEC-RAS models, is very efficient in terms of preprocessing (introducing boundary conditions in numerical models) and computing time. Our analyses show, however, that in the case of small and deep rivers with a rectangular cross-sectional shape, linear formula, Equation (1), can generate large errors in calculating river seepage, consequently falsifying simulation outcomes and their interpretations.
The base of a numerical model of seepage is two-dimensional finite element mesh (composed of nodes and elements) representing the modelled region. The accuracy of the solution depends on mesh resolution. Equation (53) is solved in each mesh node as a result of a numerical solution of a system of algebraic equations containing a coefficient matrix with a wide bandwidth of non-zero elements (the form of the bandwidth depends on the numbering of nodes and element mesh structure). Next, river–aquifer water fluxes are calculated numerically based on the derivatives of the simulated nodal total head. In the AHF model, river seepage calculation does not require construction of a spatial grid. The aquifer water flux is calculated directly as a solution of Laplace equation in separated sub-flow domains based on the data set listed in Table 1. On the other hand, the analytical model was elaborated based on a number of simplifying assumptions, in relation to the geometry of the riverbed, continuity of the flow area (assumption of no flow between P.3 and P.4 as well as between P.1 and P.3), soil layers configuration (sediment, aquifer) as well as hydraulic soil properties (only isotropic). The numerical model is free of the aforementioned limitations. It should be emphasised, however, that the application of the AHF model requires calculation of seepage through the streambed (Qbott) that is approximated with a finite series, Equation (38) of length N + 1, ( Q b o t t ( N ) ) , in which (Ak), k = 1,2,…, N + 1 is calculated from a system of 2N + 2 linear algebraic equations, Equation (36). Elements of series ( λ k ) ,   k = 1 , can be calculated from non-linear algebraic Equation (22) for instance by means of the Newton method. The convergence of the finite series ( Q b o t t ( N ) ) should be verified for increasing N . Any computing environment (allowing for solving a system of linear equations and non-linear algebraic equation), for example MATLAB software, can be used to implement the AHF model. Considering the computational complexity of the analysed models, the DM model based on Equation (1) has a zero run time, and the AHF and SEEP2D models have comparable calculation times. The advantage of the AHF model is the elimination of the time-consuming construction of the computational grid, mapping the flow domain, and ensuring acceptable assessment of river–aquifer seepage.

4. Conclusions

The simulations conducted for the exemplary river–aquifer system permitted drawing the following conclusions:
Water exchange assessed with the AHF analytical model assuming a number of simplifications can be considered the first approximation of volumetric water exchange within the exemplary river–aquifer system. When compared to the Darcy-like model DM used in many hydrogeological applications—Equation (1)—it still proves to be much more accurate.
  • The AHF model is convenient because of the simple set of data needed to solve the problem and simplicity of implementation in any computing environment.
  • The AHF model errors (estimated as a difference in total flow Qtot calculated with the AHF and SEEP2D models) depend on the “depth to width” ratio of water in the riverbed, and on the exchange flow direction-drainage or infiltration to/from the riverbed. They are in a, range of 11 to 16%, and are significantly lower compared to the DM model based on Equation (1) in which the errors are in a range of 40 to 48%.
  • A limitation of the AHF model applicability is its geometry—a rectangular-shaped riverbed cross-section followed by the same shape of the sediment layer under its bottom and alongside its bank. Overestimation of Qtot (AHF) over Qtot (SEEP2D) can be explained by restrictive assumption of horizontal flow in P.3 assumed in the AHF model.
  • For small and deep rivers, neglect of flow through the banks (as in the DM model) leads to significant errors in the total flow estimate.
Further work on the development of the AHF model will be aimed at the elimination its most restrictive assumptions.

Author Contributions

Conceptualisation, M.N., G.S., M.G.-Ł, and D.M.-Ś.; Formal analysis, M.N.; Methodology, M.N.; Software, M.N. and D.M.-Ś.; Validation, G.S. and M.G.-Ł; Visualization, G.S.; Writing—original draft, M.N., G.S., M.G.-Ł, and D.M.-Ś.; and Writing—review and editing, G.S. and M.G.-Ł. All authors have read and agreed to the published version of the manuscript.

Funding

This research was funded by Narodowe Centrum Nauki (The National Science Centre, Poland), grant number OPUS 2016/21/B/ST10/03042.

Conflicts of Interest

The authors declare no conflict of interest.

References

  1. Krause, S.; Hannah, D.M.; Fleckenstein, J.H.; Heppell, C.M.; Kaeser, D.; Pickup, R.; Pinay, G.; Robertson, A.L.; Wood, P.J. Inter-disciplinary perspectives on processes in the hyporheic zone. Ecohydrology 2011, 4, 481–499. [Google Scholar] [CrossRef] [Green Version]
  2. Boano, F.; Camporeale, C.; Revelli, R.; Ridolfi, L. Sinuosity-driven hyporheic exchange in meandering rivers. Geophys. Res. Lett. 2006, 33. [Google Scholar] [CrossRef]
  3. Jekatierynczuk-Rudczyk, E. Strefa hyporeiczna, jej funkcjonowanie i znaczenie. Kosmos 2007, 56, 181–196. [Google Scholar]
  4. Zieliński, P.; Jekatierynczuk-Rudczyk, E. Dissolved organic matter transformation in the hyporheic zone of a small lowland river. Oceanol. Hydrobiol. Stud. 2010, 39, 97–103. [Google Scholar] [CrossRef]
  5. Boano, F.; Harvey, J.W.; Marion, A.; Packman, A.I.; Revelli, R.; Ridolfi, L.; Wörman, A. Hyporheic flow and transport processes: Mechanisms, models, and biogeochemical implications. Rev. Geophys. 2014, 52, 603–679. [Google Scholar] [CrossRef]
  6. Harvey, J.; Gooseff, M. River corridor science: Hydrologic exchange and ecological consequences from bedforms to basins. Water Resour. Res. 2015, 51, 6893–6922. [Google Scholar] [CrossRef] [Green Version]
  7. Ward, A.S. The evolution and state of interdisciplinary hyporheic research. Wiley Interdiscip. Rev. Water 2016, 3, 83–103. [Google Scholar] [CrossRef]
  8. Peralta-Maraver, I.; Reiss, J.; Robertson, A.L. Interplay of hydrology, community ecology and pollutant attenuation in the hyporheic zone. Sci. Total Environ. 2018, 610–611, 267–275. [Google Scholar] [CrossRef] [Green Version]
  9. Schmadel, N.M.; Ward, A.S.; Wondzell, S.M. Hydrologic controls on hyporheic exchange in a headwater mountain stream. Water Resour. Res. 2017, 53, 6260–6278. [Google Scholar] [CrossRef]
  10. Hendriks, D.M.D.; Okruszko, T.; Acreman, M.; Grygoruk, M.; Duel, H.; Buijse, T.; Schutten, J.; Mirosław-Świątek, D.; Henriksen, H.J.; Sanches-Navarro, R.; et al. Policy Discussion Paper “Bringing groundwater to the surface” | REFORM. D7.7 Policy Discuss. Pap. No. 2, Reform Proj. 2015. Available online: https://www.reformrivers.eu/policy-discussion-paper-bringing-groundwater-surface (accessed on 6 May 2020).
  11. Trauth, N.; Schmidt, C.; Maier, U.; Vieweg, M.; Fleckenstein, J.H. Coupled 3-D stream flow and hyporheic flow model under varying stream and ambient groundwater flow conditions in a pool-riffle system. Water Resour. Res. 2013, 49, 5834–5850. [Google Scholar] [CrossRef]
  12. Siergieiev, D.; Ehlert, L.; Reimann, T.; Lundberg, A.; Liedl, R. Modelling hyporheic processes for regulated rivers under transient hydrological and hydrogeological conditions. Hydrol. Earth Syst. Sci. 2015, 19, 329–340. [Google Scholar] [CrossRef] [Green Version]
  13. Brunner, P.; Therrien, R.; Renard, P.; Simmons, C.T.; Franssen, H.J.H. Advances in understanding river-groundwater interactions. Rev. Geophys. 2017, 55, 818–854. [Google Scholar] [CrossRef]
  14. Sophocleous, M. Interactions between groundwater and surface water: The state of the science. Hydrogeol. J. 2002, 10, 52–67. [Google Scholar] [CrossRef]
  15. Grodzka-Łukaszewska, M.; Nawalany, M.; Zijl, W. A Velocity-Oriented Approach for Modflow. Transp. Porous Media 2017, 119, 373–390. [Google Scholar] [CrossRef] [Green Version]
  16. Shen, X.; Lampert, D.; Ogle, S.; Reible, D. A software tool for simulating contaminant transport and remedial effectiveness in sediment environments. Environ. Model. Softw. 2018, 109, 104–113. [Google Scholar] [CrossRef]
  17. Rushton, K.R.; Tomlinson, L.M. Possible mechanisms for leakage between aquifers and rivers. J. Hydrol. 1979, 40, 49–65. [Google Scholar] [CrossRef]
  18. Rushton, K.R.; Rathod, K.S. Horizontal and vertical components of flow deduced from groundwater heads. J. Hydrol. 1985, 79, 261–278. [Google Scholar] [CrossRef]
  19. Nawalany, M.; Grodzka-Łukaszewska, M.; Sinicyn, G. Consistency of water balances of surface waters and groundwater (in Polish). Monogr. Water Manag. Comm. Pol. Acad. Sci. 2014, 1, 237–246. [Google Scholar]
  20. Carslaw, H.S.; Jaeger, J.C. Conduction of Heat in Solids (Paperback); Oxford University Press: Oxford, UK, 1959; pp. 353–387. [Google Scholar]
  21. Harr, M.E. Groundwater and Seepage; American Association for the Advancement of Science (AAAS): Washington, DC, USA, 1963. [Google Scholar] [CrossRef]
  22. Polubarinova-Koch, P.I. Theory of Ground Water Movement; Princeton University Press: Princeton, NJ, USA, 1962. [Google Scholar] [CrossRef]
  23. Polubarinova-Koch, P.I. Theory of the Motion of Ground Water (in Russian); Nauka: Moscow, Russian, 1977. [Google Scholar]
  24. Paniconi, C.; Putti, M. Physically based modeling in catchment hydrology at 50: Survey and outlook. Water Resour. Res. 2015, 51, 7090–7129. [Google Scholar] [CrossRef] [Green Version]
  25. Cousquer, Y.; Pryet, A.; Flipo, N.; Delbart, C.; Dupuy, A. Estimating River Conductance from Prior Information to Improve Surface-Subsurface Model Calibration. Groundwater 2017, 55, 408–418. [Google Scholar] [CrossRef]
  26. Antontsev, S.N.; Díaz, J.I. On the coupling between channel level and surface ground-water flows. Pure Appl. Geophys. 2008, 165, 1511–1530. [Google Scholar] [CrossRef]
  27. Nawalany, M. Mathematical Modeling of River-Aquifer Interactions, Report SR; HR Wallingford: Wallingford, UK, 1993. [Google Scholar]
  28. Nawalany, M. Two- and three-dimensional models of groundwater flow. Phys. Appl. 1999. [Google Scholar]
  29. Saeedpanah, I.; Golmohamadi Azar, R. New Analytical Expressions for Two-Dimensional Aquifer Adjoining with Streams of Varying Water Level. Water Resour. Manag. 2017, 31, 403–424. [Google Scholar] [CrossRef]
  30. Nawalany, M. The Study of Groundwater with the Use of Dynamic System Theory; Publishing Office of Warsaw University of Technology: Warsaw, Poland, 1984. [Google Scholar]
  31. Ebel, B.A.; Mirus, B.B.; Heppner, C.S.; VanderKwaak, J.E.; Loague, K. First-order exchange coefficient coupling for simulating surface water-groundwater interactions: Parameter sensitivity and consistency with a physics-based approach. Hydrol. Process. 2009, 23, 1949–1959. [Google Scholar] [CrossRef]
  32. Bansal, R.K.; Lande, C.K.; Warke, A. Unsteady groundwater flow over sloping beds: Analytical quantification of stream-aquifer interaction in presence of thin vertical clogging layer. J. Hydrol. Eng. 2016, 21. [Google Scholar] [CrossRef]
  33. Jones, N.L. SEEP2D Primer, Environmental Modeling Research Laboratory; Brigham Young University: Provo, UT, USA, 1999. [Google Scholar]
  34. Van Genuchten, M.T. A Closed-form Equation for Predicting the Hydraulic Conductivity of Unsaturated Soils. Soil Sci. Soc. Am. J. 1980, 44, 892–898. [Google Scholar] [CrossRef] [Green Version]
Figure 1. Cross section through different shapes of riverbed: (a) Lbott >> Lbank; and (b) Lbott ≈ Lbank.
Figure 1. Cross section through different shapes of riverbed: (a) Lbott >> Lbank; and (b) Lbott ≈ Lbank.
Water 12 01792 g001
Figure 2. River–aquifer system of the Upper Biebrza River in Poland and outline of the exemplary system geometry used in the AHF and SEEP2D models.
Figure 2. River–aquifer system of the Upper Biebrza River in Poland and outline of the exemplary system geometry used in the AHF and SEEP2D models.
Water 12 01792 g002
Figure 3. Flow subdomains: P.1 and P.2, in aquifer under river sediments; P.3, river sediments right/left of the river bank; P.4, river sediments under river bottom; (y, z), general system of coordinates.
Figure 3. Flow subdomains: P.1 and P.2, in aquifer under river sediments; P.3, river sediments right/left of the river bank; P.4, river sediments under river bottom; (y, z), general system of coordinates.
Water 12 01792 g003
Figure 4. Convergence of finite series Q b o t t ( N ) for six scenarios (Φ*Hr).
Figure 4. Convergence of finite series Q b o t t ( N ) for six scenarios (Φ*Hr).
Water 12 01792 g004
Figure 5. Specific discharge fields and hydraulic head distribution in the river–aquifer system calculated for H r = 28.5 m and Φ * = 27.0 m with (a) AHF and (b) SEEP2D models (both with the assumption of no flow between subregions P.3–P.4 and P.1–P.3).
Figure 5. Specific discharge fields and hydraulic head distribution in the river–aquifer system calculated for H r = 28.5 m and Φ * = 27.0 m with (a) AHF and (b) SEEP2D models (both with the assumption of no flow between subregions P.3–P.4 and P.1–P.3).
Water 12 01792 g005
Figure 6. River recharge/discharge by groundwater seepage for different scenarios and models. Negative values of the river seepage indicate that the river is of water-gaining type; positive values occur when the river is of water-losing-type.
Figure 6. River recharge/discharge by groundwater seepage for different scenarios and models. Negative values of the river seepage indicate that the river is of water-gaining type; positive values occur when the river is of water-losing-type.
Water 12 01792 g006
Figure 7. Errors of river recharge calculated with the AHF analytical model and estimated with the Darcy’s law model (DM) related to the values calculated by the SEEP2D model.
Figure 7. Errors of river recharge calculated with the AHF analytical model and estimated with the Darcy’s law model (DM) related to the values calculated by the SEEP2D model.
Water 12 01792 g007
Figure 8. River bottom seepage for different scenarios and models.
Figure 8. River bottom seepage for different scenarios and models.
Water 12 01792 g008
Figure 9. SEEP2D calculated free groundwater table, piezometric head contours, and flow lines for (a) Hr = 28,5 m and Φ* = 27 m and (b) Hr = 25,5 m and Φ* = 27 m.
Figure 9. SEEP2D calculated free groundwater table, piezometric head contours, and flow lines for (a) Hr = 28,5 m and Φ* = 27 m and (b) Hr = 25,5 m and Φ* = 27 m.
Water 12 01792 g009
Figure 10. Two sided Qbank seepage contribution to total river flow Qtot (calculated with the SEEP2D and AHF models) as a function of the ratio of water depth (Hr–Da–ds) to river width (2Wr).
Figure 10. Two sided Qbank seepage contribution to total river flow Qtot (calculated with the SEEP2D and AHF models) as a function of the ratio of water depth (Hr–Da–ds) to river width (2Wr).
Water 12 01792 g010
Table 1. Parameters of river sediments, riverbed, and riparian aquifer in the exemplary system.
Table 1. Parameters of river sediments, riverbed, and riparian aquifer in the exemplary system.
VariableValues
Da20 m
ka0.000116 m/s
Wrs16 m
Wr4 m
ds5 m
ks0.00001 m/s
Hr25.5–28.5 m
Φ*27 m
Table 2. Impact of finite element mesh size on calculated river recharge/discharge through groundwater seepage (A, 1185 nodes; B, 4564 nodes; C, 522 nodes).
Table 2. Impact of finite element mesh size on calculated river recharge/discharge through groundwater seepage (A, 1185 nodes; B, 4564 nodes; C, 522 nodes).
Scenario No.Qtot [×10−6 m3/s/m]B/AC/A
MESH
ABC
1−39.7−39.5−39.80.99511.0027
2−27.4−27.3−27.60.99551.0057
3−14.2−14.1−14.20.99561.0026
414.914.815.00.99491.0026
530.330.230.50.99671.0050
646.446.246.60.99611.0041
Table 3. One side Q b a n k ,    Q b o t t , and Q t o t assessed by means of the analytical (AHF) and numerical (SEEP2D) models (with assumption of no flow between P.3–P.4 and P.1–P.3).
Table 3. One side Q b a n k ,    Q b o t t , and Q t o t assessed by means of the analytical (AHF) and numerical (SEEP2D) models (with assumption of no flow between P.3–P.4 and P.1–P.3).
Scenario No.Φ*[m]Hr[m]Qbank, Qbott, Qtot [×10−6 m3/s/m] 1
AHFSEEP2D
QbankQbottQtotQbankQbottQtot
127.025.5−6.89−10.76−17.65−5.39−10.56−15.95
227.026.0−4.80−7.17−11.97−3.94−7.04−10.98
327.026.5−2.51−3.58−6.09−2.13−3.52−5.65
427.027.52.713.586.292.403.525.92
527.028.05.647.1712.815.027.0412.06
627.028.58.7710.7619.537.8710.5618.43
1 negative values indicate the river is recharged by riparian aquifer groundwater through river sediments and positive values indicate the infiltration of river water into river sediments (riparian aquifer is recharged from the river).
Table 4. Calculated river recharge/discharge through groundwater seepage.
Table 4. Calculated river recharge/discharge through groundwater seepage.
Scenario No.Hr [m]River Recharge/Discharge through Groundwater Seepage Qtot [× 10−6 m3/s/m] 2
AHFSEEP2DDM
125.5−35.3−39.7−24.0
226.0−23.9−27.4−16.0
326.5−12.2−14.2−8.00
427.512.614.98.00
528.025.630.316.0
628.539.146.424.0
2 negative values indicate the river recharge by seepage from river sediments, while positive values indicate the infiltration of river water to river sediments.

Share and Cite

MDPI and ACS Style

Nawalany, M.; Sinicyn, G.; Grodzka-Łukaszewska, M.; Mirosław-Świątek, D. Groundwater–Surface Water Interaction—Analytical Approach. Water 2020, 12, 1792. https://doi.org/10.3390/w12061792

AMA Style

Nawalany M, Sinicyn G, Grodzka-Łukaszewska M, Mirosław-Świątek D. Groundwater–Surface Water Interaction—Analytical Approach. Water. 2020; 12(6):1792. https://doi.org/10.3390/w12061792

Chicago/Turabian Style

Nawalany, Marek, Grzegorz Sinicyn, Maria Grodzka-Łukaszewska, and Dorota Mirosław-Świątek. 2020. "Groundwater–Surface Water Interaction—Analytical Approach" Water 12, no. 6: 1792. https://doi.org/10.3390/w12061792

APA Style

Nawalany, M., Sinicyn, G., Grodzka-Łukaszewska, M., & Mirosław-Świątek, D. (2020). Groundwater–Surface Water Interaction—Analytical Approach. Water, 12(6), 1792. https://doi.org/10.3390/w12061792

Note that from the first issue of 2016, this journal uses article numbers instead of page numbers. See further details here.

Article Metrics

Back to TopTop