New Analytical Models for Flow Induced by Pumping in a Stream‐Aquifer System: A New Robin Boundary Condition Reflecting Joint Effect of Streambed Width and Storage

A Robin type boundary condition (BC), commonly adopted at stream‐aquifer interface, excludes a term associated with streambed accounting for the effects of streambed storage and width. This study presents two new analytical models for describing confined flow induced by pumping in a stream‐aquifer system. One model considers a single‐zone aquifer and treats the streambed as a lagging Robin BC with a time lag parameter related to the effects. The other considers a two‐zone aquifer consisting of aquifer and streambed zones. A Dirichlet BC for stream water level is specified at the edge of the streambed. The time‐domain solutions of both models are developed to describe spatiotemporal drawdown and temporal stream filtration/depletion rate (SDR). The finite element solutions (FESs) of both models are also built. Results suggest the lag time equals the half of the squared streambed width divided by the streambed hydraulic diffusivity. The effects of streambed width and storage on SDR should be considered when their lumped parameter exceeds 0.1. Neglecting their effects causes 25% difference in SDR when the lumped parameter equals 10. Based on the FESs, the use of the lagging Robin BC takes nearly a tenth of computing time of obtaining accurate steady‐state SDR for the simulation of the two‐zone aquifer. In addition, the present solutions agree to a field SDR experiment conducted by Hunt et al. (2001, https://doi.org/10.1111/j.1745‐6584.2001.tb02310.x). To conclude, this study presents two new models for describing groundwater flow in a stream‐aquifer system and explores the joint effect of streambed width and storage on SDR.


Introduction
In recent years, there is a large amount of research interest in dealing with the problems of stream-aquifer interaction (e.g., Frei & Gilfedder, 2015;Ghasemizade et al., 2019;Grinevskiy et al., 2018;Huang & Yeh, 2015;Lamontagne et al., 2015;Liang et al., 2018;Xian et al., 2019;Zipper et al., 2018Zipper et al., , 2019. Groundwater extraction from wells in an aquifer induces water filtration from an adjacent stream and may influence ecosystems due to a drop in the stream water level (Foglia et al., 2013). The stream depletion rate (SDR), defined as the ratio of a stream water filtration rate to a well discharge rate, can be used in assessing the impact of stream depletion. Huang et al. (2018) provided an extensive review on existing analytical solutions of SDR and categorized those solutions according to aquifer types (i.e., confined, unconfined, and leaky aquifers), flow dimensions (i.e., two-, three-, and quasi-three dimensions), well types (i.e., vertical, horizontal, and radial collector wells), stream representations (i.e., source terms, Dirichlet boundary, and Robin boundary), and stream channel geometry (i.e., rectilinear, strip, V-shaped, and meandering streams). Their review reveals the study for the influence of streambed storage from a single stream on SDR is unexplored using analytical approaches. Sun and Zhan (2007) studied the stream depletion induced by pumping in a confined aquifer between two parallel streams. They found that the streambed storage will affect the early stage of transient flow in the streambeds if the streambed thickness is not small as compared to the distance between two streams. Xian et al. (2020) mentioned that streambed storage is a key factor controlling the exchange flux between surface water and groundwater. They estimated a transient streambed hydraulic conductivity (K-sb) through successive inversion of flood wave responses and found the transient K-sb being significantly underestimated when neglecting the streambed storage effect. Therefore, the influence of streambed storage on the behavior of SDR may be of research interest in the development of stream-aquifer system models for water resource investigation and management.
Streambed may be clogged due to the deposition of fine particles such as silt and clay in the bed of stream. The Robin boundary condition (BC) is commonly specified at the clogged streambed between the stream and aquifer and is expressed as where s is aquifer drawdown due to pumping near a stream, x is the Cartesian coordinate normal to the stream, K is aquifer hydraulic conductivity, K ′ and w are, respectively, the streambed hydraulic conductivity and width, and the ratio K ′ /w is streambed leakance. Hantush (1965) presented an analytical solution of SDR subject to the influence of clogged streambed defined by equation 1, which in fact neglects the streambed storage. His solution reduces to the solution of Glover and Balmer (1954) when the streambed leakance approaches infinity, which leads equation 1 to the Dirichlet BC. Hunt (1999) also derived an analytical solution to SDR by treating the stream as a line source term in the flow equation. His solution is in fact exactly the same as the Hantush solution (Sun & Zhan, 2007). Sun and Zhan (2007) considered that the groundwater flow is two-dimensional (2D) in the aquifer while the flow is one-dimensional normal to the straight stream in the streambed. They investigated the water flux subject to streambed storage in a flow system of an aquifer bounded by two parallel streams. Their transient solutions however are in the Laplace domain due to the complexity of the flow system. Werner et al. (2006) investigated groundwater extraction with considering the dynamics of stream-aquifer interaction in a tropical catchment and revealed the mix of small-scale streambed and large-scale aquifer leads to much computing time due to fine grids, which may limit the use of numerical solutions when analytical solutions are not available. Thus, there may be a need to develop a new analytical model with a BC to reflect the effects of both streambed storage and permeability on the flow in a stream-aquifer system.
Recently, the lagging theory originated from the heat flow area was applied to the problems associated with pumping tests. Lin and Yeh (2017) introduced two time lag parameters in Darcy's law for describing lagging influence on flow toward a pumping well in a leaky aquifer. They indicated the time lags could be attributed to inertial effect and noninterconnected pore in the aquifer. Later, Lin et al. (2019) revealed the time lags may account for the effect of capillary suction on water table decline induced by pumping in an unconfined aquifer.
When filtration water flows from stream toward a pumping well, one might expect that the flow should take time to go through the clogged streambed. On the basis of the lagging theory, we therefore introduce a time lag parameter denoted as τ in equation 1 to handle the time lag effect induced by the streambed. The new BC is referred to as lagging Robin BC and expressed as (Lin et al., 2019;Lin & Yeh, 2017) For further exploring the physical insight of τ, two new analytical models are developed to simulate confined flow for pumping in a stream-aquifer system. One model considers a flow equation for a single-zone aquifer and regards the streambed as the lagging Robin BC. On the other hand, the other considers a two-zone aquifer consisting of aquifer and streambed zones. Each zone has its own transient 2D flow equation. In addition, the streambed zone has a Dirichlet BC specified for stream water level at the edge of the streambed. The time-domain solution of each model for the aquifer drawdown and SDR is derived by the methods of the Laplace transform and Fourier cosine transform. An analytical expression relating the τ to the hydraulic properties of the streambed is developed based on both solutions in the Laplace and Fourier domains. The finite element solutions (FESs) of the models are also built. Computation efficiency of the FESs is examined. Sensitivity analysis of SDR in response to the change in each hydraulic parameter is discussed. In addition, the present solutions are coupled with an optimization algorithm for the determination of hydraulic parameters when analyzing the SDR data measured by Hunt et al. (2001) near Doyleston Drain, New Zealand. Figure 1a shows a schematic diagram for a pumped confined aquifer near a rectilinear stream under the lagging Robin BC. The potentiometric surface before pumping is assumed at the same elevation as the stream water level. Consider both aquifer and streambed are homogeneous and isotropic. Bouwer (1969, p. 126) mentioned that the Dupuit assumption of neglecting vertical flow can be used when following two conditions are met. One is aquifer thickness beneath a partially penetrating stream less than three times of the stream width. The other is that the ratio of the distance between the well and the centerline of the stream to the stream width exceeds 10. A partially penetrating stream can be regarded as a fully penetrating stream under these two conditions. Table 1 displays the notations used in the text.

Mathematical Models
The governing equation describing 2D transient flow induced by pumping in the confined aquifer is expressed as where Q is pumping rate, T and S are, respectively, the aquifer transmissivity and storage coefficient, (x, y) is the Cartesian coordinate, and t is elapsed time since pumping starts. For a single zone aquifer, the Lagging Robin BC defined as equation 2 is specified at clogged streambed between the stream and aquifer. The initial condition is written as where the apostrophe represents streambed. The continuity requirements of drawdown and flux at the streambed-aquifer interface can be respectively denoted as The initial condition and remote boundary condition are, respectively, lim y→±∞ s ′ ¼ 0: Kollet and Zlotnik (2003, p. 97) mentioned that the pumping rate is commonly on the order of or below 0.01 m 3 /s and streamflow rate on the order of 1 m 3 /s is not uncommon. Since the effect of stream Fourier cosine parameter 10.1029/2019WR026352 Water Resources Research filtration on the streamflow is very small, thus the Dirichlet condition representing the constant stream water level is adopted and expressed as Define the following dimensionless variables and parameters: where the subscript D indicates a dimensionless symbol; S s and S ′ s are specific storages for the aquifer and streambed, respectively. With equation 11, the single-zone aquifer model, equations 2-5b, expressed in a dimensionless form is Similarly, the model for a two-zone aquifer in a dimensionless form can also be written as equations 12-14b and

Solution 1 for Single-Zone Aquifer Model
The single-zone aquifer model (i.e., equations 12-15) is solved by the methods of Fourier cosine transform and Laplace transform. The derivation is given in Appendix A. The semianalytical solution of the model is expressed as , p is the Laplace parameter, and ω is the Fourier cosine parameter. After taking the inverse Fourier cosine and Laplace transforms to equations 21a and 21b, the analytical solution of the dimensionless drawdown can be written as where F 1 , F 2 , and G are, respectively, defined as equations A14c-A14e in Appendix A. The numerical integrations for the integrals in equations 22a and 22b and other places can be done by Gaussian quadrature (Gerald & Wheatley, 2004). Each integrand has trigonometric functions and exhibits oscillation with infinite roots. Newton's method is first used to find roots. The area between each of two consecutive roots is evaluated using 16-term Gaussian quadrature formula. The integral integrated from zero to infinity is then transformed to an infinite series. The series converges very fast due to the exponential term in the integrand.
In this study, we evaluate the integrals using the Mathematica function NIntegrate (Wolfram Research Inc., 2018).
According to Darcy's law and the drawdown solution of equation 22a, the solution for SDR from the stream can be estimated by integrating flux along the boundary and dividing the result by the pumping rate and The analytical solution of SDR with the derivation given in Appendix B can be expressed as where the integrand reduces to 2+2/α for r = 0 by L'Hospital rule. Note that the SDR solution depends on both α for the streambed leakance effect and τ D for the time lag effect.

Solution 2 for Two-Zone Aquifer Model
The methods of Fourier cosine transform and Laplace transform are used to develop the solution of the two-zone aquifer model (i.e., equations 12-20 except 15). The derivation is presented in Appendix C. The solution for the drawdown in the streambed written in the Laplace and Fourier cosine domains is Similarly, the solutions for the drawdown in the aquifer on the left-hand side (LHS) and right-hand side (RHS) of the pumping well are, respectively, written as . Taking the inverse Fourier cosine and Laplace transforms to equations 25a-25c leads to the analytical solution of the dimensionless drawdown expressed as where F ′ , F 1 , F 2 , G ′ , and G are, respectively, defined as equations C10d-C10h in Appendix C.
With the derivation in Appendix D, the analytical solution of SDR defined as where ασ = κβ, σ=α ¼ w 2 D β=κ, and σ = w D β, which is a lumped parameter associated with streambed width and storage. The integrand equals 1+1/α for r = 0 by L'Hospital rule. The value of σ reflects the joint effect of the streambed width and storage on SDR as discussed in section 2.5. Note that equation 27 relies only on dimensionless parameters α and σ. Solution 1 of SDR, equation 24, depends on α and the dimensionless lag time τ D . The next section will discuss the relation of α, σ, and τ D based on Solutions 1 and 2.

Relation of Time Lag Parameter and Streamed Parameters
Let us set the second RHS term in equation 15 in the single-zone model equaling the streambed storage term in the two-zone model written as where Factor 2 is due to the symmetry of spatial drawdown distribution with respect to y D = 0. The LHS represents the total quantity of filtration water subject to the effect of time lag while the RHS accounts for the total quantity of water released from streambed storage.
Substituting equations 21a and 25a into 29 leads to , the dimensionless time lag parameter τ D , κ/w D , and lumped parameter w D β have following relation: or τ = w 2 S ′ /(2T ′ ) in a dimensional form. With equations 15 and 30, the lagging Robin BC in a dimensional form becomes Compared with equation 1, the second RHS term of equation 31 in the single-zone model indicates the time lag effect is equivalent to the joint effect of streambed width and storage in the two-zone model.

Sensitivity Analysis
The normalized sensitivity coefficient S i that evaluates change in drawdown or SDR in response to variation in each hydraulic parameter can be expressed as (Liou & Yeh, 1997) where P i is the magnitude of the ith parameter and X is either Solution 1 or Solution 2. The finite difference approach to equation 32 can be defined as with a small increment ΔP i = 10 −3 P i .

The Lagging Robin Boundary Condition and Effects of Streambed Storage and Width
This section demonstrates the applicability of the lagging Robin BC. First of all, the deviation caused by assuming streambed width as infinitesimal is examined. Consider the steady-state dimensionless drawdown s D,1 derived by the first RHS term in equation 22a for Solution 1 and equation 26b for Solution 2. The value of w D is set to unity, which represents an extreme case of close distance between stream and pumping well, to enhance the effect of the streambed width of the two-zone aquifer. Figure 2 shows spatial distributions of s D,1 at the streambed-aquifer interface (i.e., x D = 0) and dimensionless drawdown contours in the streambed and aquifer zones plotted by Solutions 1 and 2 for three scenarios of κ = 10 −3 , 10 −2 , 10 −1 regarding different ratios of streambed conductivity to aquifer one. Both solutions give close predictions of s D,1 with insignificant difference except Scenario 3 of κ = 10 −1 , which leads to minor difference in s D,1 near x D = 0 and 0 ≤ y D ≤ 10 −1 . Accordingly, we may conclude Solution 1 developed based on the lagging Robin BC gives very close prediction in steady-state aquifer drawdown to Solution 2 in spite of assuming streambed width as infinitesimal.
The following example demonstrates the behavior of SDR subject to the joint effect of streambed width and storage (i.e., wS ′ s or in dimensionless form w D β). Note that both Solutions 1 and 2 of SDR (equations 24 and 27, respectively) depend on the lumped parameters w D β (i.e., σ) and κ/w D (i.e., α) on the basis of equation 30. The SDR solution of Hunt (1999), which has the Robin BC of equation 1 and relies on only κ/w D , is chosen to make comparison with Solutions 1 and 2. The solution of Glover and Balmer (1954) with treating the stream as a Dirichlet boundary is also chosen. Cases 1-3 stand for w D β = 10, 1, and 0.1, respectively, based on a fixed ratio κ/w D = 1. Figure 3 displays temporal SDR distributions predicted by those four solutions for these three cases and drawdown contours for Cases 1 and 2 predicted by Solutions 1 and 2 with the default values w D = 0.1, t D = 1, and κ = 0.1. Both Solutions 1 and 2 agree well in the drawdown contours demonstrated in Panel (a) for Case 1 and Panel (b) for Case 2. The figure demonstrates the SDR increases with time and finally approaches unity (i.e., SDR = 1), indicating the stream filtration rate equals the pumping rate. The predicted SDR curves by both solutions have large time delays in Case 1 when t D ≤ 100 (e.g., t ≤ 6 hr), small time lag in Case 2 when t D ≤ 20, and almost no time lag in Case 3 as compared with Hunt (1999) solution when d = 50 m and T/S = 10 6 m 2 /day. The difference in the SDRs is because Solution 1 developed based on the lagging Robin BC treats the streambed width as infinitesimal while Solution 2 for the two-zone aquifer considers the streambed as a zone of finite width. The stream filtration water predicted by Solution 2 with a streambed zone therefore takes more time to flow toward the pumping well than that by Solution 1. On the other hand, Solution 1 yields a larger SDR than Solution 2 at the same pumping time. This is because Solution 2 has an additional water from the streambed storage contributed to the well pumping in addition to the water supply from the stream and aquifer storage provided by both Solutions 1 and 2.
It is worth noting that the SDR predicted by Solution 2 in Case 1 is close to zero at t D = 1 (t = 3.6 min) shown in Figure 3 although the drawdown cone already reaches the streambed-aquifer interface exhibited in Panel (a), indicating the water released from the streambed storage and aquifer storage meets the pumping demand. On the other hand, more water filtrates from the stream in Case 2 for Solution 2 and less water yielded from the streambed storage and aquifer storage, resulting in a smaller drawdown cone shown in Panel (b) than that in Panel (a). The Hunt solution neglecting streambed storage effect in Case 1 leads to the largest difference of 32% in the predicted SDR at t D = 2 compared with Solution 2 and 25% with Solution 1. For Case 3, the Hunt solution and Solutions 1 and 2 give almost the same SDR prediction because both time lag effect in Solution 1 and streambed storage effect in Solution 2 are very small and negligible.

Sensitivity Analysis
This section explores the influence of the change in each parameter on SDR behavior predicted by Solutions 1 and 2. Figure 4 shows normalized sensitivity coefficients S i for SDR to each of parameters T, S, K ′ , w, τ in Panel (a) for Solution 1 and T, S, K ′ , w, S ′ s in Panel (b) for Solution 2 when the pumping well is close to the stream (w/d = 1). The default values used in this study are w = 5 m, b = 20 m, T = 4 m 2 /hr, S = 10 −4 , K ′ = 2 × 10 −2 m/hr, S ′ s ¼ 5 × 10 −5 m −1 , and τ= 1.88 min obtained by τ ¼ w 2 S ′ s = 2K ′ À Á based on equation 30. Figure 5 displays the S i for SDR predicted by Solutions 1 and 2 but the pumping well located at some distance from the stream (w/d = 0.1). These two figures exhibit a negative S i to S ′ s or τ, indicating SDR during early pumping period decreases with increasing S ′ s because of streambed storage effect for Solution 2 and time lag effect for Solution 1. The K ′ and w are key parameters in affecting SDR. The curves of K ′ and w are symmetrical with slight difference in |S i |, indicating they are highly correlated. It is worth noting the opposite trends of S i to T in Figures 4 and 5. A larger T leads to earlier starting time of SDR > 0 because the drawdown cone takes less time to reach the stream, leading to a positive S i to T before t = 7 × 10 −3 hr in Figure 4 and during the entire period in Figure 5. The negative S i to T after t = 7 × 10 −3 hr in Figure 4 results from close distance between the well and stream (w/d = 1). A larger T causes more water from aquifer storage to the pumping well and thus decreases SDR. These results indicate the streambed storage effect or lagging effect decreases for a large distance between the well and stream, which accords with the conclusion made in Section 2.5 that both effects are negligible when w D β ≤ 0.1.

Application of Solutions 1 and 2 to Field Experiment
Hunt et al. (2001) conducted a field experiment for measuring stream filtration due to pumping in a confined aquifer adjacent to Doyleston Drain near 40 km south of Christchurch in New Zealand. The average aquifer thickness is 20 m. The formation consists of unconsolidated sand and gravel. The distance between the drain and pumping well is 55 m. The pumping was operated with a constant rate of 63 m 3 /hr for 10 hr. The streamflow was measured by two weirs, one of which is located 200 m upstream from the pumping well and the other is located 200 m downstream. This 400 m reach contains most of the filtration across the drain. The field SDR equals the ratio of the difference in the measured streamflow rates to the pumping rate as illustrated in Figure 6. One can refer to Hunt et al. (2001) for detailed description. The drain is a shallow stream of 2.5 m width overlying the aquifer. The effect of stream width W on SDR is minor when d/W ≥ 1 (d/W = 22 in this case) (Butler et al., 2001;Hunt, 2008). The drain can therefore be considered as a rectilinear boundary.
Simulation of SDR relies only on parameters T, S, K ′ /w, τ for Solution 1, T, S, K ′ /w, wS ′ s for Solution 2, and T, S, K ′ /w for the Hunt (1999) solution according to the dimensional forms of these three solutions. The values of K ′ /w and wS ′ s can therefore be estimated, but the streambed width w cannot. Note that Hunt et al. (2001) did not provide any information about w. Hunt et al. (2001) determined the aquifer parameters of T = 57.24 m 2 /hr and S = 2.11 × 10 −3 using the Hunt solution and early drawdown data before the boundary effect. They suggested determining the streambed properties by curve fitting of the SDR data because the response of drawdown is relatively insensitive to the change in streambed properties of K ′ /w and wS ′ s . With the estimates of T and S, Table 2 displays K ′ /w and τ estimated by Solution 1, K ′ /w and wS ′ s by Solution 2, and K ′ /w by the Hunt solution via coupling the Levenberg-Marquardt algorithm in the Mathematica function FindFit. The Hunt solution is developed based on the Robin BC, which does not consider the streambed storage effect. Figure 6a shows the temporal SDR distributions predicted by these three solutions with the associated parameters. Solution 1 gives the best fit with the smallest standard error of where e j represents each of ten differences in predicted and measured SDRs and v, the number of the parameters, is 3 for the Hunt solution and 4 for Solutions 1 and 2. The Hunt solution gives a larger SEE than the others because of excluding the streambed storage influence ( Table 2). Note that the estimated τ in this experiment setting may reflect the joint effects of the streambed storage and the aquifer storage beneath the drain. The latter effect should be considered because the aquifer thickness under the drain exceeds three times of the drain width (i.e., 20 m >3 × 2.5 m) (Bouwer, 1969, p. 126). This example justifies the use of the lagging Robin BC to SDR problem of field stream-aquifer system.

Computation Efficiency of Finite Element Solution
The finite element method is employed to solve the single-zone aquifer and two-zone aquifer models. The former considers the aquifer zone in the range 0 ≤ x D ≤ 3 and −4 ≤ y D ≤ 4 with the lagging Robin boundary at x D = 0 and the no-flow boundary at the other sides. The latter has the same aquifer zone as the former but the streambed zone is in

Water Resources Research
The value of SDR approaches unity when the flow reaches steady state. This condition can be used to examine the accuracy of the FES when using different size of triangle elements. Figure 7 displays the steady-state SDRs predicted by the FES when the total number of nodal points (i.e., NP) increases. The single-zone aquifer model gives good steady-state SDR of 0.99 even for the coarse grids of NP = 2,209 in Figure 7a. The two-zone aquifer model, however, yields prediction for the SDR larger than unity for the coarse grid of NP ≤ 2,290. The estimated SDR fluctuates with increasing NP, indicating the need of more nodal points allocated within and near the streambed zone. The FES yields the SDR of 0.99 when NP = 14,655 using very fine grid in the streambed zone and coarse grid in the aquifer zone shown in Figure 7b.
Define the computing time ratio as TR = CT 2 /CT 1 where CT is the computing time of the FES and subscripts 1 and 2 represent the single-zone and two-zone models, respectively. Figure 7 shows that the TR is close to 10 for the same predicted SDR from those two models, indicating the two-zone aquifer model takes nearly one order larger computing time than the other model. The CT 2 increases with more nodal points required for longer reach of stream filtration to discretize the streambed zone of the two-zone aquifer. The present single-zone model with the lagging Robin BC can be an efficient alternative to the two-zone aquifer model in predicting SDR.

Concluding Remarks
This study develops two new analytical models for confined flow induced by pumping in a stream-aquifer system. One considers a single-zone aquifer with a clogged streambed specified as a new Robin boundary condition by introducing a time lag parameter τ to reflect the joint effect of streambed width and storage. The other model considers a two-zone aquifer consisting of the aquifer and streambed zones with a constant stream water level specified as the Dirichlet condition of zero drawdown. Each zone has its own transient 2D flow equation. The analytical solutions (i.e., Solutions 1 and 2) of the models are developed for describing spatiotemporal drawdown distribution and temporal SDR distribution. The FESs of the models are also constructed to examine the accuracy of SDR prediction and required computation time. Sensitivity analysis of the response of SDR to the variation in each hydraulic parameter is investigated. On the basis of the present solutions, new findings are described below.
With the relation of τ = w 2 S ′ /(2T ′ ) (or τ D ¼ w 2 D β= 2κ ð Þ in dimensionless form), both Solutions 1 and 2 yield close predictions of drawdown and SDR with minor difference. The magnitude of w D β in Solutions 1 and 2 reflects the joint effect of streambed width and storage on SDR. When w D β > 0.1, the effect is significant. Solution 1 with the lagging Robin BC yields a smaller SDR than the one with the non-lagging Robin BC. When w D β < 0.1, the effect is negligible. The lagging Robin BC reduces to equation 1, the Robin BC. The sensitivity analysis reveals the streambed storage effect is significant when the pumping well is close to the stream. In addition, Solutions 1 and 2 are coupled with the Levenberg-Marquardt algorithm and applied to the field stream-aquifer system near Doyleston Drain, New Zealand. The lagging Robin BC based on infinitesimal width leads to better fit to the field SDR data than the treatment of regarding the streambed as a finite-width zone. These results lead to the conclusion that the single-zone model with the lagging Robin BC can be applied to estimate the SDR and reflect the joint effect of streambed width and storage.
The steady-state SDRs calculated by the FESs for the present two models have implications in the development of finite difference/element solutions. The use of the lagging Robin BC in the single-zone model allows coarse grids in the aquifer zone for the FES and thus saves computing time and yields good prediction in steady-state SDR. On the other hand, insufficient nodal points in discretizing the streambed zone in the two-zone model causes inaccurate prediction of the steady-state SDR. The lagging Robin BC is indeed a very good substitute for the streambed zone in developing finite difference/element solutions. parameter p. The transformed model becomes an ordinary differential equation (ODE) with two boundary conditions expressed as can be divided into two ODEs with two continuity requirements of drawdown and flux, written as Note that equation A12b is derived by integrating equation A7, expressed as ∫ A general formula for the inverse Laplace transform to a multiple-value function of p with a branch cut in the complex plane can be expressed as Resj p¼p n ; (A13a) is the imaginary unit, F is either b s D;1 or b s D;2 (i.e., equations 21a or 21b), f is the time-domain result for F, p 0 is a branch point, p n is the nth pole, p + = p 0 +r 2 (cosπ+isinπ) and p − = p 0 +r 2 (cosπ − isinπ) are integral paths, respectively, right above and below the branch cut located at the real axis from p = p 0 to p = − ∞, Resj p¼p n is the residue for the nth pole, and N is the number of poles. For this case we have the branch point p 0 = − ω 2 , one pole p 1 = 0, and N = 1. Note that p + and p − lead to the same result of −ω 2 − r 2 but different expressions of λ = ir for p = p + and λ = − ir for p = p − based on λ ¼ ffiffiffiffiffiffiffiffiffiffiffiffiffi ffi p þ ω 2 p . With these relations, we can have the inverse Laplace transforms to b s D;1 and b s D;2 , written as Gdr for x D ≥ 1; (A14b) where λ 1 = r 2 +ω 2 and λ 2 = τ D ω 2 − 1. Note that the results in equations A14c and A14d for ω = 0 as well as A14e for ω = 0 and r = 0 are obtained by L'Hospital rule. The formulas for the inverse Fourier cosine transform to equations A14a and A14b are expressed as The Laplace transform convertse s ′ D x D ; ω; t D ð Þintob s ′ D x D ; ω; p ð Þ , and ∂e s ′ D =∂t D in equation C1 into pb s ′ D . The transformed model can therefore be expressed as equations A8 and A10-A12b with replacingb s D byb s D;2 in A8, and Solving the model yields the Laplace-domain solution, expressed as equations 25a-25c.