Groundwater–Surface Water Interaction—Analytical Approach

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.


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): where W r is the width of the riverbed (m); k s is the hydraulic conductivity of river sediments (m/s); d s is the thickness of river sediments under the river bottom (m); Φ * is the piezometric head in the aquifer next to the riverbed (m); and H r 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 modelsthe 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.

The Flow System
Water exchange between the river and aquifer occurs both through the river bottom and river banks. For rivers where L bott >> L bank (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: D a is the thickness of the aquifer beneath river sediments (m), W rs is the half-width of river sediments (m), d s is the thickness of river sediments below the river bottom (m), W r is the half-width of the riverbed (m), L ar is the one-side length of the riparian aquifer (m), D ar 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), h s 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 bott is the vertical seepage through the river bottom per 1 meter of the river stretch segment from/to flow domain P.2 (m 3 /s/m), Q bank is the horizontal seepage through the river bank per 1 meter of the river stretch segment from/to flow domain P.3 (m 3 /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.

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: 1.
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.

2.
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.

3.
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. 4.
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. 5.
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: with boundary conditions: left boundary is the horizontal component of specific discharge at the vertical boundary between P.1 and P.2 and piezometric head continuous, i.e., right boundary is the imposed piezometric head (1st type b.c.), i.e., lower boundary is the impervious boundary (2nd type b.c.), i.e., upper boundary is the no flow boundary (2nd type b.c.), i.e., 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 where 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 rs , z) = Φ * and, consistently, ∂ Φ a (y=W rs ,z) i.e., ∂ Φ a (y, z=D a ) ∂z = 0, µ i must be equal to Notice also that total outflow from P.1 into the riparian aquifer is equal to 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: with boundary conditions: left boundary is the no flow boundary (2nd type b.c.), i.e., at the vertical boundary between P.2 and P.1, and also piezometric head is continuous, i.e., lower boundary is the impervious boundary (2nd type b.c.), i.e., upper boundary is the seepage through a layer of river sediments of P.4 to river bottom, i.e., 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 whereas its derivatives are equal to Parameters Substituting Equation (17) and Equation (19) to Equation (20) provides . . can be calculated from the non-linear algebraic equation where 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 bott can be calculated from Equation (19) by integrating specific discharge q z over top horizontal boundary of subregion P.2, i.e., 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 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) 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 rs − W r ) and added to Equation (26), resulting in Ultimately, values of (2N + 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 = mD a N , m = 0, 1 . . . , N for z into Equations (26) and (27) leads to a set of (2N + 2) linear algebraic equations, where infinite series are approximated with their finite counterparts Hyperbolic sine and cosine in Equations (30) and (31) need some pre-calculation to avoid rising exponents, i.e., sinh α = exp α [1 − exp (−2α)]/2 and cosh α = exp α [1 + exp (−2α)]/2. By denotinĝ Equations (30) and (31) can be written shortly as or in the block-matrix form for as where ] is the right hand side vector.
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 right boundary is the constant piezometric head (1st type b.c.), i.e., lower boundary is the (assumed) no flow boundary between P.1 and P.3 (2nd type b.c.), i.e., 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., 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 bank , 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 where Φ(y) = 1 h(y) h(y) 0 Φ s (y, z) dz 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 The general solution to this equation can be proposed in a simple form Substitution of parabola, Equation (46) to Equation (45) results in A = 0, and therefore 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: 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 + For , and therefore one-side For Φ * > H r , Q bank < 0, which indicates inflow to the river from P.
After moving back to the original global of coordinates (y o = 0, z o = 0), Equations (49) and (50) can be shortly written as: 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:

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 yy ∂h ∂y 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.

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 H r 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 H r 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 H r 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 (H r < 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 (H r > Φ * ), 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 Q bank , Q bott , and Q tot 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). Table 3. One side Q bank , Q bott , and Q tot assessed by means of the analytical (AHF) and numerical (SEEP2D) models (with assumption of no flow between P.3-P. 4  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).
Flow across the river bank (Q bank ) 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 (Q bott ) 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)  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 bott , 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. Q bott differences between the two models depend on the number N and range from 0.04 × 10 −6 m 3 /s/m (for N = 400, scenarios 3 and 4) to 0.2 × 10 −6 m 3 /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 tot , 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.    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-H r , Φ * , but also permits calculation of both components of the specific discharge (q y and q z ) 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 Q bott value decreases by 7% from 29.5 × 10 −6 m 3 /s/m for a drainage river to 27.5 × 10 −6 m 3 /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, Q bank accounts for approximately 40% of total flow ( Figure 10). Figure 10. Two sided Q bank seepage contribution to total river flow Q tot (calculated with the SEEP2D and AHF models) as a function of the ratio of water depth (H r -D a -d s ) to river width (2W r ). 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 (Q bott ) that is approximated with a finite series, Equation (38)

Conflicts of Interest:
The authors declare no conflict of interest.