Analysis of groundwater flow and stream depletion in L-shaped fluvial aquifers

Understanding the head distribution in aquifers is crucial for the evaluation of groundwater resources. This article develops a model for describing flow induced by pumping in an L-shaped fluvial aquifer bounded by impermeable bedrocks and two nearly fully penetrating streams. A similar scenario for numerical studies was reported in Kihm et al. (2007). The water level of the streams is assumed to be linearly varying with distance. The aquifer is divided into two subregions and the continuity conditions of the hydraulic head and flux are imposed at the interface of the subregions. The steady-state solution describing the head distribution for the model without pumping is first developed by the method of separation of variables. The transient solution for the head distribution induced by pumping is then derived based on the steady-state solution as initial condition and the methods of finite Fourier transform and Laplace transform. Moreover, the solution for stream depletion rate (SDR) from each of the two streams is also developed based on the head solution and Darcy's law. Both head and SDR solutions in the real time domain are obtained by a numerical inversion scheme called the Stehfest algorithm. The software MODFLOW is chosen to compare with the proposed head solution for the L-shaped aquifer. The steady-state and transient head distributions within the L-shaped aquifer predicted by the present solution are compared with the numerical simulations and measurement data presented in Kihm et al. (2007).

1 Introduction

Groundwater is an important water resource for agricultural, municipal and industrial uses. The planning and management of water resources through the investigation of groundwater flow is one of the major tasks for practicing engineers. The aquifer type and shape are important factors influencing the groundwater flow. Many studies have been devoted to the development of analytical models for describing flow in finite aquifers with a rectangular boundary (e.g., Chan et al., 1976, 1977; Daly and Morel-Seytoux, 1981; Latinopoulos, 1982, 1984, 1985; Corapcioglu et al., 1983; Lu et al., 2015), a wedge-shaped boundary (Chan et al., 1978; Falade, 1982; Holzbecher, 2005; Yeh et al., 2008; Chen et al., 2009; Samani and Zarei-Doudeji, 2012; Samani and Sedghi, 2015; Kacimov et al., 2016), a triangle boundary (Asadi-Aghbolaghi et al., 2010), a trapezoidal-shaped boundary (Mahdavi and Seyyedian, 2014), or a meniscus-shaped domain (Kacimov et al., 2017). So far, the case of re-entrant angle (L-shaped) boundaries has been treated analytically in different fields such as torsion of elastic bars (Kantorovich and Krylov, 1958), head fluctuation problems for tidal aquifers (Sun, 1997; Li and Jiao, 2002), and heat conduction in plates (Mackowski, 2011). However, none of the cited papers deals with pumping or stream depletion problems.

Many studies focused on the development of numerical approaches to model groundwater flow in an aquifer with irregular domain and various types of boundary conditions. The rapid increase in the computing power of PCs enables the numerical models to handle the groundwater-flow problems with complicated geometric shapes and/or a heterogeneous aquifer. Numerical methods such as finite element methods (FEMs) and finite difference methods (FDMs) are very commonly used in engineering simulations or analyses. For the application of FEMs, Taigbenu (2003) solved the transient flow problems based on the Green element method for multi-aquifer systems with arbitrary geometries. Kihm et al. (2007) used a general multidimensional hydrogeomechanical Galerkin FEM to analyze three-dimensional (3-D) problems of saturated–unsaturated flow and land displacement induced by pumping in a fluvial aquifer in the Yongpoong 2 Agriculture District, Gyeonggi-do, Korea. The domain of the aquifer is L-shaped and bounded by streams and impermeable bedrock. Their mathematical model was developed by considering the unsaturated flow and solid skeleton deformation, and thus a set of four coupled nonlinear equations in terms of one head variable and three displacement variables was derived. The Galerkin FEM was employed to simulate the steady-state spatial distribution of the hydraulic head before aquifer pumping and then the distributions of the hydraulic head and land displacement vector after 1-year pumping. Their simulation results were compared and validated with the field measurements of the hydraulic head and vertical displacement in the transient case. However, solving the simultaneous nonlinear equation required a lot of work which was difficult to apply because of the extensive calculation time and resource requirements for long-term simulations.

The FDMs have been widely utilized in the groundwater problems too. Mohanty et al. (2013) evaluated the performances of the finite difference groundwater model MODFLOW and the computational model artificial neural network (ANN) in the simulation of groundwater level in an alluvial aquifer system. They compared the results with field-observed data and found that the numerical model is suitable for long-term predictions, whereas the ANN model is appropriate for short-term applications. Serrano (2013) illustrated the use of Adomian's decomposition method to solve a regional groundwater-flow problem in an unconfined aquifer bounded by the main stream on one side, two tributaries on two sides, and an impervious boundary on the other side. He demonstrated an application to an aquifer bounded by four streams with a deep excavation inside where the head was kept constant. Jafari et al. (2016) incorporated Terzaghi's theory of one-dimensional (1-D) consolidation with MODFLOW to evaluate groundwater flow and land subsidence due to heavy pumping in a basin aquifer in Iran. So far, many computer codes developed based on either FDMs (e.g., FTWORK and MODFLOW), FEMs (e.g., AQUIFEM-N, BEMLAP, FEMWATER, and SUTRA), or boundary element methods (e.g., BEMLAP) have been employed to simulate a variety of groundwater-flow problems (Loudyi et al., 2007).

On the other hand, analytical solutions are convenient and powerful tools to explore the physical insight of groundwater-flow systems. The head solution is capable of predicting the spatiotemporal distribution of the drawdown at any location within the simulation time and the stream depletion rate (SDR) solution can estimate the stream filtration rate at any instance at a specific location in the groundwater-flow system. Thus, the development of analytical models for describing the groundwater flow in a heterogeneous aquifer with irregular outer boundaries and subject to various types of boundary conditions is of practical use from an engineering viewpoint. Kuo et al. (1994) applied the image well theory and Theis' equation to estimate transient drawdown in an aquifer with irregularly shaped boundaries. The aquifer is an oil reservoir bounded by three tortuous faults. However, the number of the image wells should be largely increased if the aquifer boundary is asymmetric and rather irregular. Insufficient number of the image wells might result in poor results or even divergence (Matthews et al., 1954). Read and Volker (1993) presented analytical solutions for steady seepage through hillsides with arbitrary flow boundaries. They used the least squares method to estimate the coefficients in a series expansion of the Laplace equation. Li et al. (1996) extended the results of Read and Volker (1993) in solving the two-dimensional (2-D) groundwater flow in porous media governed by Laplace's equation involving arbitrary boundary conditions. The solution procedure was obtained by means of an infinite series of orthonormal functions. Additionally, they also introduced a method, called the image-recharge method, to establish the recurrence relationship of the series coefficients. Patel and Serrano (2011) solved nonlinear boundary value problems of multidimensional equations by Adomian's method of decomposition for groundwater flow in irregularly shaped aquifer domains. Mahdavi and Seyyedian (2014) developed a semi-analytical solution for hydraulic head distribution in trapezoidal-shaped aquifers in response to diffusive recharge of constant rate. The aquifer was surrounded by four fully penetrating and constant-head streams. Kacimov et al. (2016) used the Strack–Chernyshov model to investigate the unconfined groundwater flows in wedge-shaped promontories with accretion along the water table and outflow from a groundwater mound into draining rays. Huang et al. (2016) presented 3-D analytical solutions for hydraulic head distributions and SDRs induced by a radial collector well in a rectangular confined or unconfined aquifer bounded by two parallel streams and no-flow boundaries. Currently, the distribution of groundwater-flow velocity in a circular meniscus aquifer was investigated analytically by the theory of holomorphic functions and numerically by FEM (Kacimov et al., 2016).

Groundwater pumping near a stream in a fluvial aquifer may cause the dispute of stream water rights, impact the aquatic ecosystem in the stream, and result in water allocation or management problems for agriculture, industry, and municipality. The impacts of groundwater extraction by wells should therefore be thoroughly investigated before pumping. Theoretically, numerical methods, such as FEM are good tools to simulate groundwater flow in irregularly shaped aquifer such as the approximately L-shaped aquifer in Kihm et al. (2007). However, they are generally time consuming and computationally demanding (Younger, 2007), and it is necessary to regenerate the mesh of problem domain for different sizes of L-shaped aquifers. The analytical solution can be easily applied for different sizes of L-shaped aquifers with similar boundaries and properties by replacing the length or width of the solution. Thus, an analytical solution is proposed in this study to evaluate the spatiotemporal distribution of the drawdown at any location in the L-shaped aquifer. This paper develops a 2-D mathematical model for describing the groundwater flow in an approximately L-shaped fluvial aquifer which is very close to the case of numerical simulations reported in Kihm et al. (2007). The aquifer is divided into two rectangular subregions. The aquifer in each subregion is assumed to be homogeneous but anisotropic in the horizontal plane with principal direction aligned with the borderline of the rectangular subregions. Three types of boundary conditions including constant head, linearly varying head, and no flow are adopted to reflect the physical reality at the outer boundaries of the problem domain. A steady-state solution is first developed to represent the hydraulic head distribution within the aquifer before pumping. The transient head solution of the model is then obtained using the Fourier finite sine and cosine transforms and the Laplace transform. The Stehfest algorithm is then taken to invert the Laplace-domain solution for the time-domain results. The software MODFLOW for the simulation of the 3-D groundwater flow is used to evaluate the present head solutions. The SDR solution is also derived based on the head solution and Darcy's law and then used to evaluate the contribution of filtration water from each of two streams toward the pumping well.

Figure 1 Location of the fluvial aquifer. Note that this figure is modified from Google Earth.

Figure 2 The L-Shaped fluvial aquifer with two subregions.

2 Methodology

Figure 1 shows a fluvial plain located in the Yongpoong 2 Agriculture District, Gyeonggi-do, Korea whose characteristics are reported in Kihm et al. (2007). The west side of the plain is a mountainous area, where impermeable bedrock outcrops, and the Poonggye stream flows along the east side from the southwest corner toward the northeast corner. A tributary of Poonggye stream, entering the stream with nearly a right angle, is on the north side of the plain. The Poonggye stream and its tributary are perennial streams and almost fully penetrate the fluvial aquifer system (Kihm et al., 2007). The width of Poonggye stream is about 15 m as reported in River Information Management GIS (RIMGIS, 2013).

2.1 Conceptual model

The aquifer in the district is formed by fluvial deposits with a total thickness of 6 m and consists of a clay loam aquitard with a thickness of about 2.5 m underlain by a loamy sand layer with a thickness of about 3.5 m (Kihm et al., 2007). In order to develop an analytical model for solving the groundwater flow, the domain of the aquifer in this study is approximated to be L-shaped, as delineated in Fig. 2. Notice that in Fig. 1 the solid line denotes the outer boundary of the L-shaped aquifer in this study while the dashed line represents the simulation area in the work of Kihm et al. (2007). The origin of the coordinates in Fig. 2 is at the lower left corner of point A, which is at the intersection of boundary AB (i.e., a part of Poonggye stream) and boundary AG. The boundaries of the aquifer domain along EF and FG are impermeable bedrocks and thus regarded as impermeable boundaries. The annual average heads above the bottom of the aquifer are respectively identified as 5.18, 4.06, and 5.29 m at points A, B, and D (Kihm et al., 2007). The hydraulic heads along AG and DE are assumed equal to their average head values done by Kihm et al. (2007). In other words, the boundaries along AG and DE are assumed under the constant-head condition in our mathematical model. Physically, they do not coincide with streams and therefore do not contribute to SDR as calculated in Sect. 2.5 (“Stream depletion rate”). The boundaries AB and BD are designated to represent the Poonggye stream and its tributary, respectively. Kihm et al. (2007) fixed the hydraulic heads of Poonggye stream and its tributary at annual average water stages in their numerical simulations. Thus, this study considers that the stream has a perfect hydraulic connection with the aquifer and the stream stage varies linearly with distance. The average stream flow rate of the Poonggye stream with its tributary is about 100 m³ s⁻¹ as reported by RIMGIS (2013, p. 90). Pumping wells in the conceptual model are assumed to fully penetrate the aquifer near the perennial stream AB as mentioned in Kihm et al. (2007), and therefore the hydraulic gradient in vertical direction is neglected. Todd and Mays (2005, p. 232) noticed that the discharge rates in a shallow well may range up to 500 m³ day⁻¹ (0.01 m³ s⁻¹) and the suction lifts should be less than 7 m for efficient and continuous service. Hence, the effect of pumping in a shallow well on the water table of the nearby stream is generally negligible. The annual average depth from the ground surface to the water table is 1.26 m with a spatial variation from 0.57 to 1.95 m in accordance with the average water stages in the streams AB and BD (Kihm et al., 2007). The depth for the aquifer system before pumping was estimated under the hydrostatic equilibrium condition and considering the effect of net annual average rainfall.

2.2 Mathematical model

As shown in Fig. 2, the aquifer is divided into two subregions named regions 1 and 2, and variables ${\mathit{\varphi }}_{\mathrm{1}}(x,y,t)$ and ${\mathit{\varphi }}_{\mathrm{2}}(x,y,t)$ are their corresponding hydraulic heads. Groundwater has been pumped by a total of M pumping wells in region 1 and N pumping wells in region 2, and all the pumping wells fully penetrate the aquifer. The coordinates of the kth well in region 1 and the lth well in region 2 are denoted as (x_1k,y_1k) and (x_2l,y_2l), respectively, and the pumping rate per unit thickness at the kth well is represented by Q_1k [L² T⁻¹] and that at the lth well is Q_2l [L² T⁻¹]. The governing equations describing 2-D hydraulic head distributions in region 1 and region 2 are respectively expressed as

$$\begin{array}{ll}\text{(1)}& {\displaystyle }& {\displaystyle }{K}_{x\mathrm{1}}{\displaystyle \frac{{\partial }^{\mathrm{2}}{\mathit{\varphi }}_{\mathrm{1}}}{\partial x^{\mathrm{2}}}}+K_{y\mathrm{1}}{\displaystyle \frac{{\partial }^{\mathrm{2}}{\mathit{\varphi }}_{\mathrm{1}}}{\partial y^{\mathrm{2}}}}={\displaystyle }& {\displaystyle }{S}_{s\mathrm{1}}{\displaystyle \frac{\partial {\mathit{\varphi }}_{\mathrm{1}}}{\partial t}}-\sum _{k=\mathrm{1}}^M{Q}_{\mathrm{1}k}\mathit{\delta }(x-x_{\mathrm{1}k})\mathit{\delta }(y-y_{\mathrm{1}k})\\ {\displaystyle }& {\displaystyle }\mathrm{0}\le x\le l_{\mathrm{1}},\mathrm{0}\le y\le d_{\mathrm{1}},\\ \text{(2)}& {\displaystyle }& {\displaystyle }{K}_{x\mathrm{2}}{\displaystyle \frac{{\partial }^{\mathrm{2}}{\mathit{\varphi }}_{\mathrm{2}}}{\partial x^{\mathrm{2}}}}+K_{y\mathrm{2}}{\displaystyle \frac{{\partial }^{\mathrm{2}}{\mathit{\varphi }}_{\mathrm{2}}}{\partial y^{\mathrm{2}}}}={\displaystyle }& {\displaystyle }{S}_{s\mathrm{2}}{\displaystyle \frac{\partial {\mathit{\varphi }}_{\mathrm{2}}}{\partial t}}-\sum _{l=\mathrm{1}}^N{Q}_{\mathrm{2}l}\mathit{\delta }(x-x_{\mathrm{2}l})\mathit{\delta }(y-y_{\mathrm{2}l})\\ {\displaystyle }& {\displaystyle }{l}_{\mathrm{2}}\le x\le l_{\mathrm{1}},d_{\mathrm{1}}\le y\le d_{\mathrm{2}},\end{array}$$

where K_x [L T⁻¹] and K_y [L T⁻¹] are respectively the hydraulic conductivities in the x and y direction, and S_s [L⁻¹] is the specific storage. The symbol δ represents the one-dimensional Dirac delta function [T⁻¹]. The boundary conditions for region 1 are expressed as

$$\begin{array}{}\text{(3)}& {\displaystyle }& {\displaystyle }{\mathit{\varphi }}_{\mathrm{1}}\left(\mathrm{0},y\right)=h_{\mathrm{1}}\text{for AG},\text{(4)}& {\displaystyle }& {\displaystyle }{\mathit{\varphi }}_{\mathrm{1}}\left(l_{\mathrm{1}},y\right)=h_{\mathrm{3}}+{\displaystyle \frac{h_{\mathrm{2}}-h_{\mathrm{3}}}{b_{\mathrm{2}}}}y\text{for BC},\text{(5)}& {\displaystyle }& {\displaystyle }{\mathit{\varphi }}_{\mathrm{1}}\left(x,\mathrm{0}\right)=h_{\mathrm{1}}+{\displaystyle \frac{h_{\mathrm{3}}-h_{\mathrm{1}}}{l_{\mathrm{1}}}}x\text{for AB},\text{(6)}& {\displaystyle }& {\displaystyle \frac{\partial {\mathit{\varphi }}_{\mathrm{1}}}{\partial y}}\left(x,d_{\mathrm{1}}\right)=\mathrm{0}\text{for FG}.\end{array}$$

Similarly, the boundary conditions for flow in region 2 are

$$\begin{array}{}\text{(7)}& {\displaystyle }& {\displaystyle \frac{\partial {\mathit{\varphi }}_{\mathrm{2}}}{\partial x}}\left(l_{\mathrm{2}},y\right)=\mathrm{0}\text{for EF},\text{(8)}& {\displaystyle }& {\displaystyle }{\mathit{\varphi }}_{\mathrm{2}}\left(l_{\mathrm{1}},y\right)=h_{\mathrm{3}}+{\displaystyle \frac{h_{\mathrm{2}}-h_{\mathrm{3}}}{b_{\mathrm{2}}}}y\text{for CD},\text{(9)}& {\displaystyle }& {\displaystyle }{\mathit{\varphi }}_{\mathrm{2}}\left(x,d_{\mathrm{2}}\right)=h_{\mathrm{2}}\text{for DE}.\end{array}$$

The continuity requirements of the hydraulic head and flux along the interface CF are respectively

$$\begin{array}{}\text{(10)}& {\mathit{\varphi }}_{\mathrm{1}}\left(x,d_{\mathrm{1}}\right)={\mathit{\varphi }}_{\mathrm{2}}(x,d_{\mathrm{1}})\end{array}$$

and

$$\begin{array}{}\text{(11)}& {\left.K_{y\mathrm{1}}{\displaystyle \frac{\partial {\mathit{\varphi }}_{\mathrm{1}}}{\partial y}}\right|}_{y=d_{\mathrm{1}}}={\left.K_{y\mathrm{2}}{\displaystyle \frac{\partial {\mathit{\varphi }}_{\mathrm{2}}}{\partial y}}\right|}_{y=d_{\mathrm{1}}}.\end{array}$$

In order to express the solution in dimensionless form, the following dimensionless variables or parameters are introduced:

• ${\mathit{\varphi }}_{\mathrm{1}}^{*}=({\mathit{\varphi }}_{\mathrm{1}}-h_{\mathrm{1}})/h_{\mathrm{1}}$,

• ${\mathit{\varphi }}_{\mathrm{2}}^{*}=({\mathit{\varphi }}_{\mathrm{2}}-h_{\mathrm{1}})/h_{\mathrm{2}}$,

• $t^{*}=K_{y\mathrm{1}}t/S_{s\mathrm{1}}{d}_{\mathrm{2}}^{\mathrm{2}}$,

• ${\mathit{\kappa }}_{\mathrm{1}}=K_{x\mathrm{1}}/K_{y\mathrm{1}}$,

• ${\mathit{\kappa }}_{\mathrm{2}}=K_{x\mathrm{2}}/K_{y\mathrm{2}}$,

• $x^{*}=x/l_{\mathrm{1}}$,

• $y^{*}=y/d_{\mathrm{2}}$,

• $d_{\mathrm{1}}^{\ast }=d_{\mathrm{1}}/d_{\mathrm{2}}$,

• $l_{\mathrm{2}}^{*}=l_{\mathrm{2}}/l_{\mathrm{1}}$,

• $Q_{\mathrm{1}k}^{\ast }=d_{\mathrm{2}}^{\mathrm{2}}{Q}_{\mathrm{1}k}/K_{y\mathrm{1}}{h}_{\mathrm{1}}$, and

• $Q_{\mathrm{2}l}^{\ast }=d_{\mathrm{2}}^{\mathrm{2}}{Q}_{\mathrm{2}l}/K_{y\mathrm{2}}{h}_{\mathrm{2}}$,

where ${\mathit{\varphi }}_{\mathrm{1}}^{*}$ and ${\mathit{\varphi }}_{\mathrm{2}}^{*}$ stand for the dimensionless hydraulic heads in regions 1 and 2, respectively; t^* refers to the dimensionless time during the test; κ₁ and κ₂ represent the anisotropic ratio of hydraulic conductivity in regions 1 and 2, respectively; x^* and y^*denote the dimensionless coordinates.

2.3 Steady-state solution for hydraulic head distribution

In order to compare the steady-state simulations of Kihm et al. (2007) without pumping, the steady-state solution for the hydraulic head distribution in the L-shaped aquifer is developed. Detailed derivation for the analytical solutions in steady state for regions 1 and 2 is given in Appendix A, and the results are expressed respectively as (Chu et al., 2012)

$$\begin{array}{ll}\text{(12)}& {\displaystyle }& {\displaystyle }{\mathit{\varphi }}_{\mathrm{1}}^{*}\left(x^{*},y^{*}\right)={\displaystyle }& {\displaystyle }\sum _{m=\mathrm{1}}^{\mathrm{\infty }}{\mathrm{\Delta }}_{\mathrm{1}}\left[C_{\mathrm{1}m}{E}_{\mathrm{1}}\left(m,y^{*}\right)+F_{\mathrm{1}}\left(m,y^{*}\right)\right]\sin \left({\mathit{\lambda }}_m{x}^{*}\right)\end{array}$$

and

$$\begin{array}{ll}\text{(13)}& {\displaystyle }& {\displaystyle }{\mathit{\varphi }}_{\mathrm{2}}^{*}\left(x^{*},y^{*}\right)={\displaystyle }& {\displaystyle }\sum _{n=\mathrm{1}}^{\mathrm{\infty }}{\mathrm{\Delta }}_{\mathrm{2}}\left[D_{\mathrm{2}n}{E}_{\mathrm{2}}\left(n,y^{*}\right)+F_{\mathrm{2}}\left(n,y^{*}\right)\right]\cos \left[{\mathit{\alpha }}_n\right(x^{*}-l_{\mathrm{2}}^{*}\left)\right],\end{array}$$

with

$$\begin{array}{}\text{(14)}& {\displaystyle }& {\displaystyle }{E}_{\mathrm{1}}\left(m,y^{*}\right)={\displaystyle \frac{e^{{\mathrm{\Omega }}_{\mathrm{1}m}{y}^{*}}-e^{-{\mathrm{\Omega }}_{\mathrm{1}m}{y}^{*}}}{e^{{\mathrm{\Omega }}_{\mathrm{1}m}{b}_{\mathrm{1}}^{*}}-e^{-{\mathrm{\Omega }}_{\mathrm{1}m}{b}_{\mathrm{1}}^{*}}}},\text{(15)}& {\displaystyle }& {\displaystyle }{F}_{\mathrm{1}}\left(m,y^{*}\right)={\displaystyle \frac{\mathrm{1}}{{\mathit{\lambda }}_m}}{\left(-\mathrm{1}\right)}^{m+\mathrm{1}}\left(h_{\mathrm{31}}^{*}+h_{\mathrm{23}}^{*}{y}^{*}\right),\text{(16)}& {\displaystyle }& {\displaystyle }{E}_{\mathrm{2}}\left(n,y^{*}\right)=e^{-{\mathrm{\Omega }}_{\mathrm{2}n}{y}^{*}}-e^{{\mathrm{\Omega }}_{\mathrm{2}n}(y^{*}-\mathrm{2})},\text{(17)}& {\displaystyle }& {\displaystyle }{F}_{\mathrm{2}}\left(n,y^{*}\right)={\displaystyle \frac{(-\mathrm{1}{)}^{n-\mathrm{1}}}{{\mathit{\alpha }}_n}}\left(H_{\mathrm{31}}^{*}+H_{\mathrm{23}}^{*}{y}^{*}\right),\end{array}$$

where the symbols and dimensionless variables Δ₁, Δ₂, λ_m, α_n, Ω_1m, Ω_2n, $h_{\mathrm{31}}^{\ast }$, $h_{\mathrm{23}}^{*}$, $H_{\mathrm{23}}^{*}$, and $H_{\mathrm{31}}^{*}$ are defined in Table 1. The coefficients C_1m and D_2n can be determined simultaneously by the continuity conditions of the hydraulic head and flux along the interface CF. The results are denoted as follows:

$$\begin{array}{ll}\text{(18)}& {\displaystyle }{C}_{\mathrm{1}m}=& {\displaystyle \frac{{\mathrm{\Delta }}_{\mathrm{2}}}{\mathrm{2}}}{\displaystyle \frac{K_{y\mathrm{2}}}{K_{y\mathrm{1}}}}{\displaystyle \frac{h_{\mathrm{2}}}{h_{\mathrm{1}}}}\sum _{n=\mathrm{1}}^{\mathrm{\infty }}\left[D_{\mathrm{2}n}{\left.{\displaystyle \frac{E_{\mathrm{2}}^{\prime }\left(n,y^{*}\right)}{E_{\mathrm{1}}^{\prime }\left(m,y^{*}\right)}}\right|}_{y^{*}=b_{\mathrm{1}}^{*}}\right.{\displaystyle }& {\displaystyle }\left.+{\left.{\displaystyle \frac{F_{\mathrm{2}}^{\prime }\left(n,y^{*}\right)}{E_{\mathrm{1}}^{\prime }\left(m,y^{*}\right)}}\right|}_{y^{*}=b_{\mathrm{1}}^{*}}\right]G_{\mathrm{1}}(m,n)-{\left.{\displaystyle \frac{F_{\mathrm{1}}^{\prime }\left(m,y^{*}\right)}{E_{\mathrm{1}}^{\prime }\left(m,y^{*}\right)}}\right|}_{y^{*}=b_{\mathrm{1}}^{*}}\end{array}$$

and

$$\begin{array}{ll}\text{(19)}& {\displaystyle }{D}_{\mathrm{2}n}=& {\displaystyle \frac{\mathrm{1}}{{\mathrm{\Delta }}_{\mathrm{2}}}}{\displaystyle \frac{h_{\mathrm{1}}}{h_{\mathrm{2}}}}\sum _{m=\mathrm{0}}^{\mathrm{\infty }}\left[C_{\mathrm{1}m}{\displaystyle \frac{E_{\mathrm{1}}(m,d_{\mathrm{1}}^{*})}{E_{\mathrm{2}}(n,d_{\mathrm{1}}^{*})}}+{\displaystyle \frac{F_{\mathrm{1}}(m,d_{\mathrm{1}}^{*})}{E_{\mathrm{2}}(n,d_{\mathrm{1}}^{*})}}\right]{\displaystyle }& {\displaystyle }{G}_{\mathrm{2}}\left(m,n\right)-{\displaystyle \frac{F_{\mathrm{2}}(n,d_{\mathrm{1}}^{*})}{E_{\mathrm{2}}(n,d_{\mathrm{1}}^{*})}},\end{array}$$

with

$$\begin{array}{}\text{(20)}& {\displaystyle }& {\displaystyle }{E}_{\mathrm{1}}^{\prime }\left(m,y^{*}\right)={\displaystyle \frac{\partial E_{\mathrm{1}}\left(m,y^{*}\right)}{\partial y^{*}}},\text{(21)}& {\displaystyle }& {\displaystyle }{E}_{\mathrm{2}}^{\prime }\left(n,y^{*}\right)={\displaystyle \frac{\partial E_{\mathrm{2}}\left(n,y^{*}\right)}{\partial y^{*}}},\end{array}$$

$$\begin{array}{}\text{(22)}& {\displaystyle }& {\displaystyle }{F}_{\mathrm{1}}^{\prime }\left(m,y^{*}\right)={\displaystyle \frac{\partial F_{\mathrm{1}}\left(m,y^{*}\right)}{\partial y^{*}}},\text{(23)}& {\displaystyle }& {\displaystyle }{F}_{\mathrm{2}}^{\prime }\left(n,y^{*}\right)={\displaystyle \frac{\partial F_{\mathrm{2}}\left(n,y^{*}\right)}{\partial y^{*}}},\text{(24)}& {\displaystyle }& {\displaystyle }{G}_{\mathrm{1}}(m,n)={\displaystyle \frac{{\int }_{l_{\mathrm{2}}}^{*\mathrm{1}}\sin \left({\mathit{\lambda }}_m{x}^{*}\right)\cos \left[{\mathit{\alpha }}_n\right(x^{*}-l_{\mathrm{2}}^{*}\left)\right]\text{d}x}{{\int }_{\mathrm{0}}^{\mathrm{1}}{\sin }^{\mathrm{2}}\left({\mathit{\lambda }}_m{x}^{*}\right)\text{d}x}},\text{(25)}& {\displaystyle }& {\displaystyle }{G}_{\mathrm{2}}(m,n)={\displaystyle \frac{{\int }_{l_{\mathrm{2}}}^{*\mathrm{1}}\sin \left({\mathit{\lambda }}_m{x}^{*}\right)\cos \left[{\mathit{\alpha }}_n\right(x^{*}-l_{\mathrm{2}}^{*}\left)\right]\text{d}x}{{\int }_{l_{\mathrm{2}}}^{*\mathrm{1}}{\cos }^{\mathrm{2}}\left[{\mathit{\alpha }}_n\right(x^{*}-l_{\mathrm{2}}^{*}\left)\right]\text{d}x}}.\end{array}$$

Table 1 Notations used in the text.

2.4 Transient solution for hydraulic head distribution

The semi-analytical solution of the model for transient hydraulic head distribution with the previous steady-state solution as the initial condition is developed via the methods of finite sine transform, finite cosine transform, and Laplace transform. The detailed derivation for the transient solution is given in Appendix B and the results of the dimensionless hydraulic heads in Laplace domain for regions 1 and 2 are respectively

$$\begin{array}{ll}\text{(26)}& {\displaystyle }& {\displaystyle }{\stackrel{\mathrm{̃}}{\mathit{\varphi }}}_{\mathrm{1}}^{*}\left(x^{*},y^{*},p\right)={\displaystyle }& {\displaystyle }{\mathit{\delta }}_{\mathrm{1}}\sum _{i=\mathrm{1}}^{\mathrm{\infty }}\left[w_{\mathrm{1}i}^{*}{\text{TE}}_{\mathrm{1}}\left(i,y^{*},p\right)+T_{\mathrm{1}}\left(i,y^{*},p\right)\right.\\ {\displaystyle }& {\displaystyle }\left.+T_{\mathrm{2}}\left(i,y^{*},p\right)+{\text{SQ}}_{\mathrm{1}}\left(i,y^{*},p\right)\right]\sin \left({\mathit{\lambda }}_i{x}^{*}\right)\end{array}$$

and

$$\begin{array}{ll}\text{(27)}& {\displaystyle }& {\displaystyle }{\stackrel{\mathrm{̃}}{\mathit{\varphi }}}_{\mathrm{2}}^{*}\left(x^{*},y^{*},p\right)={\displaystyle }& {\displaystyle }{\mathit{\delta }}_{\mathrm{2}}\sum _{j=\mathrm{1}}^{\mathrm{\infty }}\left[w_{\mathrm{2}j}^{*}{\text{TE}}_{\mathrm{2}}\left(j,y^{*},p\right)+T_{\mathrm{4}}\left(j,y^{*},p\right)\right.\\ {\displaystyle }& {\displaystyle }\left.+T_{\mathrm{5}}\left(j,y^{*},p\right)+{\text{SQ}}_{\mathrm{2}}\left(j,y^{*},p\right)\right]\cos \left[{\mathit{\alpha }}_j\right(x^{*}-l_{\mathrm{2}}^{*}\left)\right],\end{array}$$

with

$$\begin{array}{}\text{(28)}& {\displaystyle }{\displaystyle }{\text{TE}}_{\mathrm{1}}\left(i,y^{*},p\right)={\displaystyle \frac{e^{{\mathit{\mu }}_i(y^{*}-d_{\mathrm{1}}^{*})}-e^{-{\mathit{\mu }}_i(y^{*}+d_{\mathrm{1}}^{*})}}{\mathrm{1}-e^{-\mathrm{2}{\mathit{\mu }}_i{d}_{\mathrm{1}}^{*}}}},\end{array}$$

$$\begin{array}{ll}\text{(29)}& {\displaystyle }& {\displaystyle }{T}_{\mathrm{1}}\left(i,y^{*},p\right)={\displaystyle }& {\displaystyle \frac{\mathrm{1}}{{\mathit{\mu }}_{i}p}}\left[{\mathit{\theta }}_{\mathrm{1}}^{\mathrm{2}}{\mathit{\lambda }}_i{\left(-\mathrm{1}\right)}^i\right]\left[h_{\mathrm{31}}^{*}{e}^{-{\mathit{\mu }}_i{y}^{*}}+\left(h_{\mathrm{23}}^{*}{y}^{*}-h_{\mathrm{31}}^{*}\right)\right]\\ {\displaystyle }& {\displaystyle }-h_{\mathrm{31}}^{*}{\displaystyle \frac{{\left(-\mathrm{1}\right)}^i}{{\mathit{\lambda }}_i}}{e}^{-{\mathit{\mu }}_i{y}^{*}},\end{array}$$

$$\begin{array}{ll}\text{(30)}& {\displaystyle }& {\displaystyle }{T}_{\mathrm{2}}\left(i,y^{*},p\right)={\displaystyle }& {\displaystyle }-C_{\mathrm{1}m}\left[{\displaystyle \frac{\mathrm{1}}{\mathrm{2}}}-{\displaystyle \frac{\sin \left(\mathrm{2}{\mathit{\lambda }}_i\right)}{\mathrm{4}{\mathit{\lambda }}_i}}\right]\left[{\displaystyle \frac{-e^{-{\mathrm{\Omega }}_{\mathrm{1}i}\left(y^{*}-d_{\mathrm{1}}^{*}\right)}+e^{{\mathrm{\Omega }}_{\mathrm{1}i}\left(y^{*}-d_{\mathrm{1}}^{*}\right)}}{\left(\mathrm{1}-e^{\mathrm{2}{\mathrm{\Omega }}_{\mathrm{1}i}{d}_{\mathrm{1}}^{*}}\right)({{\mathrm{\Omega }}_{\mathrm{1}i}}^{\mathrm{2}}-{\mathit{\mu }}_i^{\mathrm{2}})}}\right.\\ {\displaystyle }& {\displaystyle }\left.+{\displaystyle \frac{{\mathit{\delta }}_{\mathrm{1}}{\left(-\mathrm{1}\right)}^i}{{\mathit{\mu }}_i^{\mathrm{2}}}}\left(h_{\mathrm{23}}^{*}{y}^{*}+h_{\mathrm{31}}^{*}-h_{\mathrm{31}}^{*}{e}^{-{\mathit{\mu }}_i{y}^{*}}\right)\right],\end{array}$$

$$\begin{array}{ll}\text{(31)}& {\displaystyle }& {\displaystyle }{\text{SQ}}_{\mathrm{1}}\left(i,y^{*},p\right)={\displaystyle }& {\displaystyle }\left\{\begin{array}{ll}{\displaystyle \frac{\mathrm{1}}{\mathrm{2}{\mathit{\mu }}_{i}p}}\sum _{k=\mathrm{1}}^M{Q}_{\mathrm{1}k}^{*}\sin \left({\mathit{\lambda }}_i{x}_{\mathrm{1}k}^{*}\right){\displaystyle \frac{\mathrm{1}}{\mathrm{1}-e^{-{\mathit{\mu }}_i{d}_{\mathrm{1}}^{*}}}}& y^{*}>y_{\mathrm{1}k}^{*}\\ \left[\begin{array}{c}{e}^{{\mathit{\mu }}_i\left(y_{\mathrm{1}k}^{*}-\mathrm{2}{d}_{\mathrm{1}}^{*}\right)}+e^{{\mathit{\mu }}_i\left(y^{*}-y_{\mathrm{1}k}^{*}-d_{\mathrm{1}}^{*}\right)}\\ +e^{{\mathit{\mu }}_i\left(y_{\mathrm{1}k}^{*}-y^{*}\right)}-e^{{\mathit{\mu }}_i\left(y^{*}-y_{\mathrm{1}k}^{*}-\mathrm{2}{d}_{\mathrm{1}}^{*}\right)}\\ -e^{{\mathit{\mu }}_i\left(y_{\mathrm{1}k}^{*}-y^{*}-d_{\mathrm{1}}^{*}\right)}-e^{-{\mathit{\mu }}_i{y}_{\mathrm{1}k}^{*}}\end{array}\right]& \\ {\displaystyle \frac{\mathrm{1}}{\mathrm{2}{\mathit{\mu }}_{i}p}}\sum _{k=\mathrm{1}}^M{Q}_{\mathrm{1}k}^{*}\sin \left({\mathit{\lambda }}_i{x}_{\mathrm{1}k}^{*}\right){\displaystyle \frac{\mathrm{1}}{\mathrm{1}-e^{-{\mathit{\mu }}_i{d}_{\mathrm{1}}^{*}}}}& y^{*}<y_{\mathrm{1}k}^{*}\\ \left[\begin{array}{c}{e}^{{\mathit{\mu }}_i(y^{*}-y_{\mathrm{1}k}^{*})}+e^{{\mathit{\mu }}_i(y_{\mathrm{1}k}^{*}-\mathrm{2}{d}_{\mathrm{1}}^{*})}\\ -e^{{\mathit{\mu }}_i\left(y^{*}+y_{\mathrm{1}k}^{*}-\mathrm{2}{d}_{\mathrm{1}}^{*}\right)}-e^{-{\mathit{\mu }}_i{y}_{\mathrm{1}k}^{*}}\end{array}\right]& \end{array}\right.,\end{array}$$

$$\begin{array}{}\text{(32)}& {\displaystyle }{\displaystyle }{\text{TE}}_{\mathrm{2}}\left(j,y^{*},p\right)={\displaystyle \frac{e^{-{\mathit{\theta }}_j{y}^{*}}-e^{{\mathit{\theta }}_j(y^{*}-\mathrm{2})}}{e^{-{\mathit{\theta }}_j{d}_{\mathrm{1}}^{*}}-e^{{\mathit{\theta }}_j(d_{\mathrm{1}}^{*}-\mathrm{2})}}},\end{array}$$

$$\begin{array}{ll}\text{(33)}& {\displaystyle }& {\displaystyle }{T}_{\mathrm{4}}\left(j,y^{*},p\right)={\displaystyle }& {\displaystyle }\left[{\displaystyle \frac{{\mathit{\theta }}_{\mathrm{2}}^{\mathrm{2}}{\mathit{\alpha }}_j{\left(-\mathrm{1}\right)}^j}{p}}+{\displaystyle \frac{S_{s\mathrm{2}}{K}_{y\mathrm{1}}{\left(-\mathrm{1}\right)}^j}{S_{s\mathrm{1}}{K}_{y\mathrm{2}}{\mathit{\alpha }}_j}}\right]\\ {\displaystyle }& {\displaystyle }\left[{\displaystyle \frac{H_{\mathrm{31}}^{*}+H_{\mathrm{23}}^{*}{y}^{*}-H_{\mathrm{21}}^{*}{e}^{{\mathit{\theta }}_j(y^{*}-\mathrm{1})}}{{\mathit{\theta }}_j^{\mathrm{2}}}}\right]+{\displaystyle \frac{H_{\mathrm{21}}^{*}}{{\mathit{\alpha }}_j}}{\displaystyle \frac{{\left(-\mathrm{1}\right)}^j}{p}}{e}^{{\mathit{\theta }}_j(y^{*}-\mathrm{1})},\end{array}$$

$$\begin{array}{ll}\text{(34)}& {\displaystyle }& {\displaystyle }{T}_{\mathrm{5}}\left(j,y^{*},p\right)={\displaystyle }& {\displaystyle }-D_{\mathrm{2}n}{\displaystyle \frac{S_{s\mathrm{2}}{K}_{y\mathrm{1}}(\mathrm{1}-l_{\mathrm{2}}^{*})}{\mathrm{2}{S}_{\text{s1}}{K}_{y\mathrm{2}}}}\left({\displaystyle \frac{e^{{\mathrm{\Omega }}_{\mathrm{2}j}(y^{*}-\mathrm{2})}-e^{-{\mathrm{\Omega }}_{\mathrm{2}j}{y}^{*}}}{{\mathit{\theta }}_j^{\mathrm{2}}-{\mathrm{\Omega }}_{\mathrm{2}j}}}\right),\end{array}$$

$$\begin{array}{ll}\text{(35)}& {\displaystyle }& {\displaystyle }{\text{SQ}}_{\mathrm{2}}\left(j,y^{*},p\right)={\displaystyle }& {\displaystyle }\left\{\begin{array}{ll}{\displaystyle \frac{\mathrm{1}}{\mathrm{2}{\mathit{\theta }}_{j}p}}\sum _{l=\mathrm{1}}^N{Q}_{\mathrm{2}l}^{*}\cos \left({\mathit{\alpha }}_j{x}_{\mathrm{2}l}^{*}\right)& y^{*}>y_{\mathrm{2}j}^{*}\\ \left(e^{-{\mathit{\theta }}_j\left(y^{*}-y_{\mathrm{2}l}^{*}\right)}-e^{-{\mathit{\theta }}_j\left(\mathrm{2}-y^{*}-y_{\mathrm{2}l}^{*}\right)}\right)& \\ {\displaystyle \frac{\mathrm{1}}{\mathrm{2}{\mathit{\theta }}_{j}p}}\sum _{l=\mathrm{1}}^N{Q}_{\mathrm{2}l}^{*}\cos \left({\mathit{\alpha }}_j{x}_{\mathrm{2}l}^{*}\right)& y^{*}<y_{\mathrm{2}j}^{*}\\ \left(e^{{\mathit{\theta }}_j(y^{*}-y_{\mathrm{2}l}^{*})}-e^{-{\mathit{\theta }}_j(\mathrm{2}-y^{*}-y_{\mathrm{2}l}^{*})}\right)& \end{array}\right.,\end{array}$$

where p is the Laplace variable and the symbols or dimensionless parameters δ₁, δ₂, α_j, λ_i, μ_i, θ₁, θ₂, θ_j, Ω_1i, Ω_2j, and $H_{\mathrm{21}}^{*}$ are introduced in Table 1.

The coefficients in Eqs. (26) and (27) are obtained via continuity requirements for the hydraulic head and flow flux at the interface CF. They can be solved simultaneously based on the following two equations

$$\begin{array}{ll}\text{(36)}& {\displaystyle }& {\displaystyle }{w}_{\mathrm{1}i}^{*}={\displaystyle \frac{K_{y\mathrm{2}}}{K_{y\mathrm{1}}}}{\displaystyle \frac{h_{\mathrm{2}}}{h_{\mathrm{1}}}}{G}_{\mathrm{1}}\left(i,j\right){\displaystyle }& {\displaystyle }\sum _{i=\mathrm{1}}^{\mathrm{\infty }}{\left.{\displaystyle \frac{{\displaystyle \begin{array}{c}{w}_{\mathrm{2}j}^{*}{\text{TE}}_{\mathrm{2}}^{\prime }\left(j,y^{*},p\right)+T_{\mathrm{4}}^{\prime }\left(j,y^{*},p\right)\\ +T_{\mathrm{5}}^{\prime }\left(j,y^{*},p\right)+{\text{SQ}}_{\mathrm{2}}^{\prime }\left(j,y^{*},p\right)\end{array}}}{{\text{TE}}_{\mathrm{1}}^{\prime }\left(j,y^{*},p\right)}}\right|}_{y^{*}=d_{\mathrm{1}}^{*}}\\ {\displaystyle }+& {\displaystyle }\sum _{i=\mathrm{1}}^{\mathrm{\infty }}{\left.{\displaystyle \frac{{\displaystyle \begin{array}{c}-T_{\mathrm{1}}^{\prime }\left(j,y^{*},p\right)+T_{\mathrm{2}}^{\prime }\left(j,y^{*},p\right)\\ +{\text{SQ}}_{\mathrm{1}}^{\prime }\left(j,y^{*},p\right)\end{array}}}{{\text{TE}}_{\mathrm{1}}^{\prime }\left(j,y^{*},p\right)}}\right|}_{y^{*}=d_{\mathrm{1}}^{*}}\end{array}$$

and

$$\begin{array}{ll}\text{(37)}& {\displaystyle }& {\displaystyle }{w}_{\mathrm{2}j}^{*}={\displaystyle \frac{h_{\mathrm{1}}}{h_{\mathrm{2}}}}{G}_{\mathrm{2}}\left(j,i\right){\displaystyle }& {\displaystyle }\sum _{j=\mathrm{1}}^{\mathrm{\infty }}{\left.{\displaystyle \frac{{\displaystyle \begin{array}{c}{w}_{\mathrm{1}i}^{*}{\text{TE}}_{\mathrm{1}}\left(i,y^{*},p\right)+T_{\mathrm{1}}\left(i,y^{*},p\right)\\ -T_{\mathrm{2}}\left(i,y^{*},p\right)+{\text{SQ}}_{\mathrm{1}}\left(i,y^{*},p\right)\end{array}}}{{\text{TE}}_{\mathrm{2}}\left(j,y^{*},p\right)}}\right|}_{y^{*}=d_{\mathrm{1}}^{*}}\\ {\displaystyle }-& {\displaystyle }\sum _{j=\mathrm{1}}^{\mathrm{\infty }}{\left.{\displaystyle \frac{{\displaystyle \begin{array}{c}{T}_{\mathrm{4}}\left(j,y^{*},p\right)+T_{\mathrm{5}}\left(j,y^{*},p\right)\\ +{\text{SQ}}_{\mathrm{2}}\left(j,y^{*},p\right)\end{array}}}{{\text{TE}}_{\mathrm{2}}\left(j,y^{*},p\right)}}\right|}_{y^{*}=d_{\mathrm{1}}^{*}},\end{array}$$

with

$$\begin{array}{}\text{(38)}& {\displaystyle }& {\displaystyle }{\text{TE}}_{\mathrm{1}}^{\prime }\left(i,y^{*},p\right)={\displaystyle \frac{\partial {\text{TE}}_{\mathrm{1}}\left(i,y^{*},p\right)}{\partial y^{*}}},\text{(39)}& {\displaystyle }& {\displaystyle }{\text{TE}}_{\mathrm{2}}^{\prime }\left(j,y^{*},p\right)={\displaystyle \frac{\partial {\text{TE}}_{\mathrm{2}}\left(j,y^{*},p\right)}{\partial y^{*}}},\text{(40)}& {\displaystyle }& {\displaystyle }{T}_{\mathrm{1}}^{\prime }\left(i,y^{*},p\right)={\displaystyle \frac{\partial T_{\mathrm{1}}\left(i,y^{*},p\right)}{\partial y^{*}}},\text{(41)}& {\displaystyle }& {\displaystyle }{T}_{\mathrm{2}}^{\prime }\left(i,y^{*},p\right)={\displaystyle \frac{\partial T_{\mathrm{2}}\left(i,y^{*},p\right)}{\partial y^{*}}},\text{(42)}& {\displaystyle }& {\displaystyle }{T}_{\mathrm{4}}^{\prime }\left(j,y^{*},p\right)={\displaystyle \frac{\partial T_{\mathrm{4}}\left(j,y^{*},p\right)}{\partial y^{*}}},\text{(43)}& {\displaystyle }& {\displaystyle }{T}_{\mathrm{5}}^{\prime }\left(j,y^{*},p\right)={\displaystyle \frac{\partial T_{\mathrm{5}}\left(j,y^{*},p\right)}{\partial y^{*}}},\text{(44)}& {\displaystyle }& {\displaystyle }S{Q}_{\mathrm{1}}^{\prime }\left(i,y^{*},p\right)={\displaystyle \frac{\partial {\text{SQ}}_{\mathrm{1}}\left(i,y^{*},p\right)}{\partial y^{*}}},\text{(45)}& {\displaystyle }& {\displaystyle }{\text{SQ}}_{\mathrm{2}}^{\prime }\left(j,y^{*},p\right)={\displaystyle \frac{\partial {\text{SQ}}_{\mathrm{2}}\left(j,y^{*},p\right)}{\partial y^{*}}}.\end{array}$$

The coefficient $w_{\mathrm{2}j}^{*}$ can be determined by substituting Eq. (36) into Eq. (37); $w_{\mathrm{1}i}^{*}$ can then be obtained once $w_{\mathrm{2}j}^{*}$ is known. The hydraulic head distributions in the real time domain can be obtained by applying a numerical Laplace inversion scheme called the Stehfest algorithm (Stehfest, 1970), to Eqs. (26) and (27).

2.5 Stream depletion rate

Pumping in an aquifer near a stream often produces water filtration from the stream toward the well (Yeh et al., 2008). Water extracted by the pumping well comes from sources such as aquifer storage and nearby streams. The extraction rate from the stream is referred to as the stream depletion rate. Since the boundaries AG and ED do not correspond to streams in the physical world and are mathematically treated as constant head because they are far from the pumping well, only the water filtration from streams AB and BD to the nearby pumping well needs to be considered. The dimensionless solutions of SDR in the Laplace domain from the stream reaches AB and BD, denoted respectively as ${\stackrel{\mathrm{̃}}{\text{SDR}}}_{\text{AB}}$ and ${\stackrel{\mathrm{̃}}{\text{SDR}}}_{\text{BD}}$, can be estimated by taking the derivatives of Eqs. (26) and (27) with respect to y and x, respectively, and then integrating along the reaches as

$$\begin{array}{}\text{(46)}& {\stackrel{\mathrm{̃}}{\text{SDR}}}_{\text{AB}}={\displaystyle \frac{q_{\text{A}}}{Q}}=-{\displaystyle \frac{\mathrm{1}}{Q}}\underset{\mathrm{0}}{\overset{l_{\mathrm{1}}}{\int }}{\left.K_{y\mathrm{1}}{\displaystyle \frac{\partial {\stackrel{\mathrm{̃}}{\mathit{\varphi }}}_{\mathrm{1}}\left(x,y,p\right)}{\partial y}}\right|}_{y=\mathrm{0}}\text{d}x\end{array}$$

and

$$\begin{array}{ll}\text{(47)}& {\displaystyle }{\stackrel{\mathrm{̃}}{\text{SDR}}}_{\mathrm{BD}}={\displaystyle \frac{q_{\text{B}}}{Q}}={\displaystyle \frac{\mathrm{1}}{Q}}(& {\displaystyle }\underset{\mathrm{0}}{\overset{d_{\mathrm{1}}}{\int }}{K}_{x\mathrm{1}}{\left.{\displaystyle \frac{\partial {\stackrel{\mathrm{̃}}{\mathit{\varphi }}}_{\mathrm{1}}\left(x,y,p\right)}{\partial x}}\right|}_{x=l_{\mathrm{1}}}\text{d}y{\displaystyle }+& {\displaystyle }\left.\underset{d_{\mathrm{1}}}{\overset{d_{\mathrm{2}}}{\int }}{K}_{x\mathrm{2}}{\left.{\displaystyle \frac{\partial {\stackrel{\mathrm{̃}}{\mathit{\varphi }}}_{\mathrm{2}}\left(x,y,p\right)}{\partial x}}\right|}_{x=l_{\mathrm{1}}}\text{d}y\right).\end{array}$$

Further, the SDR solution for streams AB and BD in the real time domain are respectively denoted as SDR_AB and SDR_BD and also obtained by the Stehfest algorithm (Stehfest, 1970). The total dimensionless stream depletion rate (SDR_T) in the time domain comes from the streams (AB and BD) and is expressed as

$$\begin{array}{}\text{(48)}& {\text{SDR}}_{\text{T}}={\text{SDR}}_{\text{AB}}+{\text{SDR}}_{\mathrm{BD}}.\end{array}$$

Since the depletion rate from constant-head boundaries AG and DE, which are far from the pumping well, can be neglected, the dimensionless storage release rate (SRR) representing the storage release rate due to compression of the aquifer matrix and expansion of groundwater in the pore space can be approximated as

$$\begin{array}{}\text{(49)}& \text{SRR}=\mathrm{1}-{\text{SDR}}_{\text{T}}.\end{array}$$

Figure 3 Contours of the hydraulic head in the L-shaped aquifer predicted by the present solution, MODFLOW, and FEM simulations with an irregular outer boundary reported in Kihm et al. (2007).

3 Comparisons of the present solution, numerical solutions, and field-observed data

3.1 Comparisons of the present solution with the MODFLOW solution

The software MODFLOW (USGS, 2005) is used to simulate the groundwater flow due to pumping in the L-shaped aquifer in the Yongpoong 2 Agriculture District with different hydraulic conductivities for the two layers. As shown in Fig. 1, region 1 has an area of 852 m × 222 m (i.e., l₁×d₁), while the area of region 2 is 297 m × 183 m (i.e., $(l_{\mathrm{1}}-l_{\mathrm{2}})\times (d_{\mathrm{2}}-d_{\mathrm{1}}))$. Thus, the total area of these two regions is 243 495 m², which is close to the area of the fluvial aquifer (246 500 m²) reported in Kihm et al. (2007). In the simulation of MODFLOW, the plane of the L-shaped aquifer is discretized with a uniform cell size of 3 m × 3 m. The aquifer thickness is 6 m and divided into two layers. The upper loam layer is 2.5 m thick and the lower sand layer is 3.5 m thick (Kihm et al. 2007). Within the aquifer domain, there is a total of 54 110 cells while the total number of cells are 42 032 and 12 078 respectively for region 1 and region 2. The boundary conditions specified for the L-shaped aquifer are the same as those defined in the mathematical model. The hydraulic heads along AG and DE are respectively h₁=5.18 m and h₂=5.29 m and the head at point B is h₃=4.06 m. Following Kihm et al. (2007), the fluvial aquifer is considered isotropic and homogeneous in the horizontal direction. In other words, the hydraulic conductivities in the x and y directions are identical in both regions 1 and 2 (i.e., $K_{x\mathrm{1}}=K_{y\mathrm{1}}=K_{x\mathrm{2}}=K_{y\mathrm{2}}=K)$. However, the aquifer is heterogeneous in the vertical direction. It has two layers with hydraulic conductivity $K_{\mathrm{U}}=\mathrm{3}\times {\mathrm{10}}^{-\mathrm{6}}$ m s⁻¹ for the upper layer and $K_{\mathrm{L}}=\mathrm{2}\times {\mathrm{10}}^{-\mathrm{4}}$ m s⁻¹ for the lower layer. The specific storage of the aquifer in both regions 1 and 2 is 10⁻⁴ m⁻¹ (Kihm et al., 2007). Consider that the pumping well P_w is located at (609 m, 9 m) in region 1 shown in Fig. 2 with a rate of 120 m³ day⁻¹ for 1-year pumping. Figure 3 shows the hydraulic head distribution obtained from MODFLOW simulations and denoted as the dotted line.

The global behavior of a multilayered aquifer may be approximated with that of an equivalent homogeneous medium, whose hydraulic conductivity in the horizontal plane K_h may be evaluated as (Charbeneau, 2000):

$$\begin{array}{}\text{(50)}& K_{\mathrm{h}}=\sum _i^m{b}_i{K}_i/\sum _i^m{b}_i,\end{array}$$

where K_i is the hydraulic conductivity in the horizontal direction for layer i, b_i is the thickness of layer i, and m is the number of the layers. Accordingly, the equivalent horizontal hydraulic conductivity K_h for the two-layered L-shaped aquifer is estimated as 1.2 × 10⁻⁴ m s⁻¹. The solid line in Fig. 3 represents the hydraulic head distribution predicted by the present solution of Eqs. (26) and (27). The head distribution predicted by the present solution agrees with that of MODFLOW simulations except in the region near the no-flow boundary FG, where it has the largest relative deviation at 2.1 %. The comparison of the head distributions indicates that the use of equivalent hydraulic conductivity in the present model is appropriate and gives fairly good predicted results.

3.2 Steady-state head distribution without pumping in the Yongpoong 2 Agriculture District and impact of aquifer anisotropy

Kihm et al. (2007) reported the steady-state hydraulic head distribution, shown in Fig. 4 by the dashed line, for the FEM simulation without groundwater pumping in the two-layered irregular aquifer. Figure 4 also shows the steady-state head distributions predicted by the present solution of Eqs. (11) and (12) denoted as the solid line and by the MODFLOW denoted as the dotted line both for the L-shaped aquifer with $K_{x\mathrm{1}}=K_{y\mathrm{1}}=K_{x\mathrm{2}}=K_{y\mathrm{2}}=\mathrm{1.2}\times {\mathrm{10}}^{-\mathrm{4}}$ m s⁻¹ (i.e., ${\mathit{\kappa }}_{\mathrm{1}}={\mathit{\kappa }}_{\mathrm{2}}=\mathrm{1})$ evaluated based on Eq. (50) and other aquifer properties mentioned in Sect. 3.1. Similar to the pumping case, the predicted head distribution of the present solution conforms to the results of the MODFLOW simulation except in the area near the no-flow boundary FG. The result of the present solution shows that the predicted contour lines of the head distribution are nearly parallel to the boundary AG and perpendicular to the boundary FG in the region x≤200 m. Moreover, the predicted heads within the regions between 500 m $\le x\le$ 852 m and 0 m $\le y\le$ 200 m are reasonably close to the FEM results, which range from 4.3 to 4.7 m as shown in Fig. 4. The groundwater flows toward point B since it has the lowest water table within the problem domain.

Figure 4 Steady-state hydraulic head contours without pumping in the Yongpoong 2 Agriculture District.

Figure 5 Steady-state hydraulic head contours in the L-shaped aquifers with three different anisotropy ratios for ${\mathit{\kappa }}_{\mathrm{1}}={\mathit{\kappa }}_{\mathrm{2}}=\mathit{\kappa }$.

Figure 5 shows the contour lines of the hydraulic head distribution for the isotropic case of ${\mathit{\kappa }}_{\mathrm{1}}={\mathit{\kappa }}_{\mathrm{2}}=\mathrm{1}$ by the solid line and for anisotropic cases of ${\mathit{\kappa }}_{\mathrm{1}}={\mathit{\kappa }}_{\mathrm{2}}=\mathrm{4}$ represented by the dashed–dotted line and ${\mathit{\kappa }}_{\mathrm{1}}={\mathit{\kappa }}_{\mathrm{2}}=\mathrm{0.25}$ by the dashed line. In these three cases, the head distributions are significantly different in the region where x≤600 m for the head ranging from 5 to 4.6 m. The largest head difference occurs near the upper boundary FG, reflecting the effects of no-flow condition and aquifer anisotropy on the flow pattern within this area.

3.3 Spatial head distributions due to pumping simulated by Kihm et al. (2007) and present solution after 1-year pumping

Figure 3 shows the spatial head distributions in the L-shaped aquifer predicted by the present solution and the MODFLOW for 1-year pumping at well P_w located at (609 m, 9 m) with a rate of 120 m³ day⁻¹. In fact, Kihm et al. (2007) reported their FEM simulations for head distributions, groundwater-flow velocity, and land displacement for 1-year pumping at the well P_w with the same pumping rate mentioned above. They referred the simulated head results as initial steady-state distributions for the case of no pumping and final steady-state distributions for the case after 1-year pumping. The aquifer configuration in their FEM simulations and the simulated head distributions denoted as a dashed line are also demonstrated in Fig. 3. The figure indicates that the present solution gives good predicted head contours near the pumping well and a reasonably good result for the head distribution in region 1 as compared to those given by Kihm et al. (2007). The head distributions predicted by the FEM solution and present solution have obvious differences in the area far away from the pumping well. Those differences may be mainly caused by the difference in the physical domain considered in the FEM solution and the simplified domain used in the present solution. In addition, the mathematical model in Kihm et al. (2007) considered the unsaturated flow and deformation of the unsaturated soil, which may also affect the head distribution after pumping. Notice that the pumping well is very close to the stream boundary AB, which is the main stream in that area and provides a large amount of filtration water to the well. Hence, both groundwater flows in region 1 for x≤300 m (near boundary AG) and in region 2 for y≥200 m (near boundaries DE and EF) are not noticeably influenced by the pumping because these two regions are far away from the pumping well.

Three piezometers, O₁, O₂, and O₃, were respectively installed at (597 m, 25 m), (594 m, 48 m), and (597 m, 204 m) as mentioned in Kihm et al. (2007) and indicated in Fig. 2. Note that O₁ was installed near the stream AB while O₃ was far away from the stream but close to the impermeable boundary FG. Figure 6 shows the temporal distributions of the hydraulic head measured at these three piezometers (i.e., H_iO, $i=\mathrm{1},\mathrm{2},\mathrm{3})$ and predicted by the FEM simulations (Kihm et al., 2007) (i.e., H_iF), present solution (i.e., H_iA), and MODFLOW simulations (i.e., H_iM). This figure indicates that the hydraulic heads predicted by the present solution has a good agreement with those simulated by Kihm et al. (2007). The largest relative differences between the temporal hydraulic head predicted by the FEM and present solution are 0.58, 0.31, and 0.51 % m for O₁, O₂, and O₃, respectively. As shown in Figs. 2 and 3, more than half of the aquifer boundary is surrounded by perennial streams and constant-head boundaries, and the pumping well is installed close to the main stream AB for stream filtration and screened in the highly permeable sand layer (Kihm et al., 2007). It is expected that the extracted water from the pumping well mainly comes from the stream AB, and the differences of the hydraulic heads from the observations (O₁, O₂, and O₃) between the present solution and FEM simulation are insignificant because the boundary geometry near these observation in both the present solution and Kihm et al. (2007) are similar. These difference may be caused by the solid skeleton deformation and unsaturated flow considered in the FEM simulation of Kihm et al. (2007). Young's modulus (E) of the sand layer is 1.9×10⁷ N m⁻² (Kihm et al., 2007), implying that the deformation of the sand layer due to pumping may be small. This is also seen in the vertical displacements reported in Kihm et al. (2007). The largest vertical displacement reported in Kihm et al. (2007) is only −0.003 m. Hence, the effect of land displacement during pumping may not significantly influence the hydraulic heads in piezometers. In addition, the unsaturated flow may slightly affect the saturated flow system. This is because the average thickness of the unsaturated zone is about 1.26 m, and it consists of a low permeability material (i.e., loam), which is 2 orders of magnitude less than that of the saturated zone (i.e., most of sand). As mentioned earlier, the pumping well is screened in the sand layer and the stream AB is the major source of extracted water. Therefore, the influence of the unsaturated zone on the saturated flow may be small. Under such hydrogeological conditions, the present solution yields a similar prediction for the temporal hydraulic head distribution as compared with those of FEM. Compared with the field observation, the differences of the predicted hydraulic head among FEM, the present solution, and MODFLOW are all less than 0.08 m at these three piezometers during 0.1 to 10 days. In addition, the largest relative differences between measured heads and predicted heads by the present solution at O₁ to O₃ during 0.1 to 5 days are respectively 1.64, 1.74, and 0.62 %, indicating that the present solution gives good predictions in the early pumping period. Moreover, the effects of unsaturated flow and land deformation on the groundwater flow in the Yongpoong aquifer are small and may be neglected. The hydraulic head at O₁ drops more than those at O₂ and O₃, whereas the former is located closer to the pumping well P_w. Because P_w is very near the stream, the extracted water will be quickly contributed from the stream and therefore the drawdown at O₁ will be soon stabilized. Figure 6 indicates that the hydraulic heads at O₁−O₃ predicted by the present solution reach steady state after t=100, 220, and 290 days, respectively.

Figure 6 Temporal distributions of the hydraulic head observed at piezometer O_i, H_iO, and H_iF are, respectively, the observations and FEM simulations reported in Kihm et al. (2007); H_iA and H_iM are the prediction of the present solution and MODFLOW, respectively, for i= 1, 2, 3.

3.4 Stream filtration in fluvial aquifer systems

Stream filtration can be considered as a problem associated with the interaction between the groundwater and surface water. The pumped water originated from the nearby stream is commonly supplied for irrigation, municipalities, and rural homes. In stream basins with several tributaries, pumping wells are often installed adjacent to the confluence of two tributaries in fluvial aquifers (Lambs, 2004).

Figure 7 Temporal distributions of SDRs, SRR, and CHRs (constant head contribution rates) due to pumping at P_w.

It is of practical interest to know the temporal SDR distributions from both streams in the Yongpoong area when subject to pumping at P_w under a rate of 120 m³ day⁻¹. The distances from P_w to the streams AB and BD are respectively 9 and 243 m. Figure 7 shows the temporal SDR distribution from each stream, indicating that SDR_AB (i.e., SDR from stream AB) is significantly larger than SDR_BD (SDR from stream BD) all the time. The steady-state values for SDR_AB and SDR_BD are respectively 0.81 and 0.19 when t≥220 days. This is due to the fact that the pumping well is closer to stream AB than stream BD and therefore water contributing to the pumping well from stream AB is much more than that from stream BD. Figure 7 also shows that the SDR_T is zero and the aquifer storage release rate SRR is unity when t≤0.01 days, indicating that the well discharges totally from the aquifer storage in the early pumping period. Once the drawdown cone reaches the stream, the SDR_T increases quickly with time while the SRR decreases continuously over the entire pumping period. It is interesting to mention that SDR_T starts to contribute more water to the pumping than SRR when t≥5 days. Finally, the SDR_T reaches unity and the SRR equals zero after t≥220 days, indicating that the aquifer system approaches steady state and all the extraction water comes from the streams.

4 Conclusions

A new semi-analytical model has been developed to analyze the 2-D hydraulic head distributions with and without pumping in a heterogeneous and anisotropic aquifer for an L-shaped domain bounded by two streams with linearly varying hydraulic heads. The method of domain decomposition is used to divide the aquifer into two regions for the development of the semi-analytical solution. The steady-state solution is first derived and used as the initial condition for the L-shaped aquifer system before pumping. The Laplace-domain solution of the model for transient head distribution in the aquifer subject to pumping is developed using the Fourier finite sine and cosine transform and the Laplace transform. The solution for SDR describing the filtration rate from two streams in an L-shaped aquifer is developed based on the head solution and Darcy's law. The Stehfest algorithm is then adopted to evaluate the time-domain results for both the head and SDR solutions in Laplace domain.

The 3-D finite difference model MODFLOW is first used to support the evaluation of the hydraulic head predictions by the present solution for the L-shaped two-layered aquifer system. The hydraulic head distributions predicted by present solutions agree fairly well over the entire aquifer except for the heads nearing the no-flow boundary. The solution for hydraulic head distribution in the L-shaped aquifer without pumping has been used to investigate the effect of anisotropic ratio (K_x∕K_y) on the steady-state flow system. It is interesting to note that the flow pattern in terms of lines of equal hydraulic head is strongly influenced by the value of anisotropic ratio for the region near the turning point of the L-shaped aquifer.

The transient solution for head distribution is employed to simulate the head distribution induced by pumping in the aquifer within the agriculture area of Gyeonggi-do, Korea. The aquifer is approximated as L-shaped in this study. The present solution delivers fairly good results in predicting the temporal hydraulic head distribution while comparing with those of FEM reported in previous study. Those simulation results seem to indicate that the effects of unsaturated flow and land displacement on the groundwater flow are not significant and may be ignorable. The largest relative differences between the measured heads and the predicted heads by the present solution at three piezometers are less than 1.74 %.

The SDR solution is first used to evaluate the steady-state SDR from each of the nearby streams for the Yongpoong aquifer subject to a specific pumping rate. The solution is also employed to determine the temporal contribution rates from the aquifer storage and the streams toward the extraction well.

Appendix A: Steady-state solution for flow in an L-shaped aquifer without pumping

On the basis of dimensionless variables and parameters defined in Sect. 2.2, Eqs. (1) and (2) can be written respectively as

$$\begin{array}{ll}\text{(A1)}& {\displaystyle }& {\displaystyle }{\mathit{\kappa }}_{\mathrm{1}}{\displaystyle \frac{b_{\mathrm{2}}^{\mathrm{2}}}{l_{\mathrm{1}}^{\mathrm{2}}}}{\displaystyle \frac{{\partial }^{\mathrm{2}}{\mathit{\varphi }}_{\mathrm{1}}^{*}}{\partial x^{*\mathrm{2}}}}+{\displaystyle \frac{{\partial }^{\mathrm{2}}{\mathit{\varphi }}_{\mathrm{1}}^{*}}{\partial y^{*\mathrm{2}}}}={\displaystyle }& {\displaystyle \frac{\partial {\mathit{\varphi }}_{\mathrm{1}}^{*}}{\partial t^{*}}}-\sum _{k=\mathrm{1}}^M{Q}_{\mathrm{1}k}^{*}\mathit{\delta }\left(x^{*}-x_{\mathrm{1}k}^{*}\right)\mathit{\delta }\left(y^{*}-y_{\mathrm{1}k}^{*}\right)\\ {\displaystyle }& {\displaystyle }\mathrm{0}\le x^{*}\le l_{\mathrm{1}}^{*},\mathrm{0}\le y^{*}\le d_{\mathrm{1}}^{*}\end{array}$$

and

$$\begin{array}{ll}\text{(A2)}& {\displaystyle }& {\displaystyle }{\mathit{\kappa }}_{\mathrm{2}}{\displaystyle \frac{d_{\mathrm{2}}^{\mathrm{2}}}{l_{\mathrm{1}}^{\mathrm{2}}}}{\displaystyle \frac{{\partial }^{\mathrm{2}}{\mathit{\varphi }}_{\mathrm{2}}^{*}}{\partial x^{*\mathrm{2}}}}+{\displaystyle \frac{{\partial }^{\mathrm{2}}{\mathit{\varphi }}_{\mathrm{2}}^{*}}{\partial y^{*\mathrm{2}}}}={\displaystyle }& {\displaystyle \frac{S_{s\mathrm{2}}}{S_{s\mathrm{1}}}}{\displaystyle \frac{K_{y\mathrm{1}}}{K_{y\mathrm{2}}}}{\displaystyle \frac{\partial {\mathit{\varphi }}_{\mathrm{2}}^{*}}{\partial t^{*}}}-\sum _{l=\mathrm{1}}^N{Q}_{\mathrm{2}l}^{*}\mathit{\delta }\left(x^{*}-x_{\mathrm{2}l}^{*}\right)\mathit{\delta }\left(y^{*}-y_{\mathrm{2}l}^{*}\right)\\ {\displaystyle }& {\displaystyle }\mathrm{0}\le x^{*}\le l_{\mathrm{2}}^{*},d_{\mathrm{1}}^{*}\le y^{*}\le \mathrm{1}.\end{array}$$

The dimensionless boundary conditions for region 1 can be expressed as

$$\begin{array}{}\text{(A3)}& {\displaystyle }& {\displaystyle }{\mathit{\varphi }}_{\mathrm{1}}^{*}\left(\mathrm{0},y^{*}\right)=\mathrm{0}\text{for boundary AG},\text{(A4)}& {\displaystyle }& {\displaystyle }{\mathit{\varphi }}_{\mathrm{1}}^{*}\left(\mathrm{1},y^{*}\right)=h_{\mathrm{31}}^{*}+h_{\mathrm{23}}^{*}{y}^{*}\text{for boundary BC},\text{(A5)}& {\displaystyle }& {\displaystyle }{\mathit{\varphi }}_{\mathrm{1}}^{*}\left(x^{*},\mathrm{0}\right)=h_{\mathrm{31}}^{*}{x}^{*}\text{for boundary AB},\text{(A6)}& {\displaystyle }& {\displaystyle \frac{\partial {\mathit{\varphi }}_{\mathrm{1}}^{*}}{\partial y^{*}}}\left(x^{*},d_{\mathrm{1}}^{*}\right)=\mathrm{0}\text{for FG},\end{array}$$

and for region 2 as

$$\begin{array}{}\text{(A7)}& {\displaystyle }& {\displaystyle \frac{\partial {\mathit{\varphi }}_{\mathrm{2}}^{*}}{\partial x^{*}}}\left(l_{\mathrm{2}}^{*},y^{*}\right)=\mathrm{0}\text{for boundary EF},\text{(A8)}& {\displaystyle }& {\displaystyle }{\mathit{\varphi }}_{\mathrm{2}}^{*}\left(\mathrm{1},y^{*}\right)=h_{\mathrm{31}}^{*}+h_{\mathrm{23}}^{*}{y}^{*}\text{for boundary CD},\text{(A9)}& {\displaystyle }& {\displaystyle }{\mathit{\varphi }}_{\mathrm{2}}^{*}\left(x^{*},\mathrm{1}\right)=h_{\mathrm{21}}^{*}\text{for boundary DE}.\end{array}$$

The continuity requirements of the hydraulic head and flux at the region interface in dimensionless form are respectively expressed as

$$\begin{array}{}\text{(A10)}& {\displaystyle }{h}_{\mathrm{1}}{\mathit{\varphi }}_{\mathrm{1}}^{*}\left(x^{*},d_{\mathrm{1}}^{*}\right)=h_{\mathrm{2}}{\mathit{\varphi }}_{\mathrm{2}}^{*}\left(x^{*},d_{\mathrm{1}}^{*}\right)\text{for segment CF}\end{array}$$

and

$$\begin{array}{ll}\text{(A11)}& {\displaystyle }& {\displaystyle }{\left.K_{y\mathrm{1}}{h}_{\mathrm{1}}{\displaystyle \frac{\partial {\mathit{\varphi }}_{\mathrm{1}}^{*}}{\partial y^{*}}}\right|}_{y^{*}=d_{\mathrm{1}}^{*}}={\displaystyle }& {\displaystyle }{\left.K_{y\mathrm{2}}{h}_{\mathrm{2}}{\displaystyle \frac{\partial {\mathit{\varphi }}_{\mathrm{2}}^{*}}{\partial y^{*}}}\right|}_{y^{*}=d_{\mathrm{1}}^{*}}\text{for segment CF}.\end{array}$$

The steady-state solution for groundwater flow in an L-shaped aquifer without pumping can be solved after removing the source/sink term in Eqs. (A1) and (A2). Multiplying Eq. (A1) by sin (λ_mx^*) and integrating it for x^* from 0 to 1 in region 1 with boundary conditions Eqs. (A3) and (A4), Eq. (A1) is then transformed to

$$\begin{array}{}\text{(A12)}& {\displaystyle }{\mathrm{\Omega }}_{\mathrm{1}m}^{\mathrm{2}}{\stackrel{\mathrm{‾}}{\mathit{\varphi }}}_{\mathrm{1}}^{*}-{\displaystyle \frac{{\partial }^{\mathrm{2}}{\stackrel{\mathrm{‾}}{\mathit{\varphi }}}_{\mathrm{1}}^{*}}{\partial y^{*\mathrm{2}}}}=-{\mathit{\kappa }}_{\mathrm{1}}{\displaystyle \frac{d_{\mathrm{2}}^{\mathrm{2}}}{l_{\mathrm{1}}^{\mathrm{2}}}}{\mathit{\lambda }}_m{\left(-\mathrm{1}\right)}^m(h_{\mathrm{31}}^{*}+h_{\mathrm{23}}^{*}{y}^{*}),\end{array}$$

with

$$\begin{array}{}\text{(A13)}& {\displaystyle }{\stackrel{\mathrm{‾}}{\mathit{\varphi }}}_{\mathrm{1}}^{*}=\underset{\mathrm{0}}{\overset{\mathrm{1}}{\int }}{\mathit{\varphi }}_{\mathrm{1}}^{*}\sin \left({\mathit{\lambda }}_m{x}^{*}\right)\text{d}{x}^{*},\end{array}$$

where ${\mathrm{\Omega }}_{\mathrm{1}m}={\mathit{\lambda }}_m\sqrt{{\mathit{\kappa }}_{\mathrm{1}}}{d}_{\mathrm{2}}/l_{\mathrm{1}}$, λ_m=mπ and $m=\mathrm{1},\mathrm{2},\mathrm{3},\mathrm{\dots }$ .

Similarly, Eq. (A2) can be transformed via multiplying Eq. (A2) by $\cos \left[{\mathit{\alpha }}_n\right(x^{*}-l_{\mathrm{2}}^{*}\left)\right]$ and integrating it for x^* from $l_{\mathrm{2}}^{*}$ to 1 in region 2 with boundary conditions Eqs. (A7) and (A8). The result is

$$\begin{array}{}\text{(A14)}& {\displaystyle }{\displaystyle }{\mathrm{\Omega }}_{\mathrm{2}n}^{\mathrm{2}}{\stackrel{\mathrm{‾}}{\mathit{\varphi }}}_{\mathrm{2}}^{*}-{\displaystyle \frac{{\partial }^{\mathrm{2}}{\stackrel{\mathrm{‾}}{\mathit{\varphi }}}_{\mathrm{2}}^{*}}{\partial y^{*\mathrm{2}}}}={\mathit{\kappa }}_{\mathrm{2}}{\displaystyle \frac{d_{\mathrm{2}}^{\mathrm{2}}}{l_{\mathrm{1}}}},\end{array}$$

with

$$\begin{array}{}\text{(A15)}& {\displaystyle }{\stackrel{\mathrm{‾}}{\mathit{\varphi }}}_{\mathrm{2}}^{*}=\underset{l_{\mathrm{2}}}{\overset{*\mathrm{1}}{\int }}{\mathit{\varphi }}_{\mathrm{2}}^{*}\cos \left[{\mathit{\alpha }}_n\left(x^{*}-l_{\mathrm{2}}^{*}\right)\right]\text{d}{x}^{*},\end{array}$$

where ${\mathrm{\Omega }}_{\mathrm{2}n}={\mathit{\alpha }}_n\sqrt{{\mathit{\kappa }}_{\mathrm{2}}}{d}_{\mathrm{2}}/l_{\mathrm{1}}$, ${\mathit{\alpha }}_n=(n-\mathrm{1}/\mathrm{2})\mathit{\pi }/(\mathrm{1}-l_{\mathrm{2}}^{*})$ and $n=\mathrm{1},\mathrm{2},\mathrm{3},\mathrm{\dots }$ .

The general solutions of Eqs. (A12) and (A14) can be written respectively as

$$\begin{array}{ll}\text{(A16)}& {\displaystyle }& {\displaystyle }{\stackrel{\mathrm{‾}}{\mathit{\varphi }}}_{\mathrm{1}}^{*}\left(m,y^{*}\right)={\displaystyle }& {\displaystyle }{C}_{\mathrm{1}m}{e}^{{\mathrm{\Omega }}_{\mathrm{1}m}{y}^{*}}+C_{\mathrm{2}m}{e}^{-{\mathrm{\Omega }}_{\mathrm{1}m}{y}^{*}}-{\displaystyle \frac{{\left(-\mathrm{1}\right)}^m}{{\mathit{\lambda }}_m}}(h_{\mathrm{31}}^{*}+h_{\mathrm{23}}^{*}{y}^{*})\end{array}$$

and

$$\begin{array}{ll}\text{(A17)}& {\displaystyle }& {\displaystyle }{\stackrel{\mathrm{‾}}{\mathit{\varphi }}}_{\mathrm{2}}^{*}\left(n,y^{*}\right)={\displaystyle }& {\displaystyle }{D}_{\mathrm{1}n}{e}^{{\mathrm{\Omega }}_{\mathrm{2}n}{y}^{*}}+D_{\mathrm{2}n}{e}^{-{\mathrm{\Omega }}_{\mathrm{2}n}{y}^{*}}-{\displaystyle \frac{{\left(-\mathrm{1}\right)}^{n-\mathrm{1}}}{{\mathit{\alpha }}_n}}(H_{\mathrm{31}}^{*}+H_{\mathrm{23}}^{*}{y}^{*}).\end{array}$$

The coefficients C_1m and C_2m in Eq. (A16) are determined by Eq. (A5) and the result is

$$\begin{array}{}\text{(A18)}& {\displaystyle }{C}_{\mathrm{2}m}=-C_{\mathrm{1}m}.\end{array}$$

Similarly, the coefficients D_1n and D_2n in Eq. (A17) are determined based on Eq. (A10) as

$$\begin{array}{}\text{(A19)}& {\displaystyle }{D}_{\mathrm{1}n}=-D_{\mathrm{2}n}{e}^{-\mathrm{2}{\mathrm{\Omega }}_{\mathrm{2}n}}.\end{array}$$

Substituting Eq. (A18) into Eq. (A16), the inversion of ${\stackrel{\mathrm{‾}}{\mathit{\varphi }}}_{\mathrm{1}}^{*}$ leads to Eq. (12) for the dimensionless hydraulic head distribution in region 1. Similarly, the inversion of ${\stackrel{\mathrm{‾}}{\mathit{\varphi }}}_{\mathrm{2}}^{*}$ for region 2 after substituting Eq. (A19) into Eq. (A17) results in Eq. (13). Based on Eqs. (A10) and (A11), the coefficients of C_1m and D_2n can be simultaneously determined and the results are respectively given in Eqs. (18) and (19).

Appendix B: Transient solutions for an L-shaped aquifer

Multiplying Eq. (A1) by sin (λ_ix^*) and integrating it for x^* from 0 to 1 in region 1 with Eqs. (A3) and (A4), Eq. (A1) can be transformed as

$$\begin{array}{ll}\text{(B1)}& {\displaystyle }& {\displaystyle }-{\mathrm{\Omega }}_{\mathrm{1}i}^{\mathrm{2}}{\stackrel{\mathrm{‾}}{\mathit{\varphi }}}_{\mathrm{1}}^{*}-{\mathit{\theta }}_{\mathrm{1}}^{\mathrm{2}}{\mathit{\lambda }}_i{\left(-\mathrm{1}\right)}^i\left(h_{\mathrm{31}}^{*}+h_{\mathrm{23}}^{*}{y}^{*}\right)+{\displaystyle \frac{{\partial }^{\mathrm{2}}{\stackrel{\mathrm{‾}}{\mathit{\varphi }}}_{\mathrm{1}}^{*}}{\partial y^{*\mathrm{2}}}}={\displaystyle }& {\displaystyle \frac{\partial {\stackrel{\mathrm{‾}}{\mathit{\varphi }}}_{\mathrm{1}}^{*}}{\partial t^{*}}}-\sum _{k=\mathrm{1}}^M{Q}_{\mathrm{1}k}^{*}\sin \left({\mathit{\lambda }}_i{x}_{\mathrm{1}k}^{*}\right)\mathit{\delta }\left(y^{*}-y_{\mathrm{1}k}^{*}\right),\end{array}$$

with

$$\begin{array}{}\text{(B2)}& {\displaystyle }{\stackrel{\mathrm{‾}}{\mathit{\varphi }}}_{\mathrm{1}}^{*}=\underset{\mathrm{0}}{\overset{\mathrm{1}}{\int }}{\mathit{\varphi }}_{\mathrm{1}}^{*}\sin \left({\mathit{\lambda }}_i{x}_{\mathrm{1}k}^{*}\right)\text{d}{x}^{*},\end{array}$$

where ${\mathrm{\Omega }}_{\mathrm{1}i}={\mathit{\lambda }}_i\sqrt{{\mathit{\kappa }}_{\mathrm{1}}}{d}_{\mathrm{2}}/l_{\mathrm{1}}$, ${\mathit{\theta }}_{\mathrm{1}}=\sqrt{{\mathit{\kappa }}_{\mathrm{1}}}{d}_{\mathrm{2}}/l_{\mathrm{1}}$, and λ_i=iπ, $i=\mathrm{1},\mathrm{2},\mathrm{3},\mathrm{\dots }$ .

Similarly, Eq. (A2) can be transformed via multiplying Eq. (A2) by cos (α_jx^*) and integrating it for x^* from $l_{\mathrm{2}}^{*}$ to 1 in region 2 with Eqs. (A7) and (A8). The result is

$$\begin{array}{ll}\text{(B3)}& {\displaystyle }& {\displaystyle }-{\mathrm{\Omega }}_{\mathrm{2}j}^{\mathrm{2}}{\stackrel{\mathrm{‾}}{\mathit{\varphi }}}_{\mathrm{2}}^{*}+{\mathit{\theta }}_{\mathrm{2}}^{\mathrm{2}}{\mathit{\alpha }}_j{\left(-\mathrm{1}\right)}^j\left(H_{\mathrm{31}}^{*}+H_{\mathrm{23}}^{*}{y}^{*}\right)+{\displaystyle \frac{{\partial }^{\mathrm{2}}{\stackrel{\mathrm{‾}}{\mathit{\varphi }}}_{\mathrm{2}}^{*}}{\partial y^{*\mathrm{2}}}}={\displaystyle }& {\displaystyle \frac{K_{y\mathrm{1}}{S}_{s\mathrm{2}}}{K_{y\mathrm{2}}{S}_{s\mathrm{1}}}}{\displaystyle \frac{\partial {\stackrel{\mathrm{‾}}{\mathit{\varphi }}}_{\mathrm{2}}^{*}}{\partial t^{*}}}-\sum _{l=\mathrm{1}}^N{Q}_{\mathrm{2}l}^{*}\cos \left({\mathit{\alpha }}_j{x}_{\mathrm{2}l}^{*}\right)\mathit{\delta }\left(y^{*}-y_{\mathrm{2}l}^{*}\right),\end{array}$$

with

$$\begin{array}{}\text{(B4)}& {\displaystyle }{\stackrel{\mathrm{‾}}{\mathit{\varphi }}}_{\mathrm{2}}^{*}=\underset{l_{\mathrm{2}}}{\overset{*\mathrm{1}}{\int }}{\mathit{\varphi }}_{\mathrm{2}}^{*}\cos \left({\mathit{\alpha }}_j{x}_{\mathrm{2}l}^{*}\right)\text{d}{x}^{*},\end{array}$$

where ${\mathrm{\Omega }}_{\mathrm{2}j}={\mathit{\alpha }}_j\sqrt{{\mathit{\kappa }}_{\mathrm{2}}}{d}_{\mathrm{2}}/l_{\mathrm{1}}$, ${\mathit{\theta }}_{\mathrm{2}}=\sqrt{{\mathit{\kappa }}_{\mathrm{2}}}{d}_{\mathrm{2}}/l_{\mathrm{1}},$ and ${\mathit{\alpha }}_j=\left(\mathrm{1}-\mathrm{1}/\mathrm{2}\right)\mathit{\pi }/(\mathrm{1}-l_{\mathrm{2}}^{*})$ for $j=\mathrm{1},\mathrm{2},\mathrm{3},\mathrm{\dots }$ . Then, taking Laplace transforms to Eq. (B1) results in

$$\begin{array}{ll}\text{(B5)}& {\displaystyle }& {\displaystyle }-{\mathrm{\Omega }}_{\mathrm{1}i}^{\mathrm{2}}{\stackrel{\mathrm{‾}}{\stackrel{\mathrm{̃}}{\mathit{\varphi }}}}_{\mathrm{1}}^{*}-{\displaystyle \frac{\mathrm{1}}{p}}{\mathit{\theta }}_{\mathrm{1}}^{\mathrm{2}}{\mathit{\lambda }}_i{\left(-\mathrm{1}\right)}^i\left(h_{\mathrm{31}}^{*}+h_{\mathrm{23}}^{*}{y}^{*}\right)+{\displaystyle \frac{{\partial }^{\mathrm{2}}{\stackrel{\mathrm{‾}}{\stackrel{\mathrm{̃}}{\mathit{\varphi }}}}_{\mathrm{1}}^{*}}{\partial y^{*\mathrm{2}}}}={\displaystyle }& {\displaystyle }p{\stackrel{\mathrm{‾}}{\stackrel{\mathrm{̃}}{\mathit{\varphi }}}}_{\mathrm{1}}^{*}-{\stackrel{\mathrm{‾}}{\mathit{\varphi }}}_{\mathrm{1}s}^{*}-{\displaystyle \frac{\mathrm{1}}{p}}\sum _{k=\mathrm{1}}^M{Q}_{\mathrm{1}k}^{*}\sin \left({\mathit{\lambda }}_i{x}_{\mathrm{1}k}^{*}\right)\mathit{\delta }\left(y^{*}-y_{\mathrm{1}k}^{*}\right),\end{array}$$

where ${\stackrel{\mathrm{‾}}{\mathit{\varphi }}}_{\mathrm{1}s}^{*}$ is the steady-state solution of region 1. Hence, Eq. (B5) can be organized as

$$\begin{array}{ll}\text{(B6)}& {\displaystyle }& {\displaystyle }-{\mathit{\mu }}_i^{\mathrm{2}}+{\stackrel{\mathrm{‾}}{\stackrel{\mathrm{̃}}{\mathit{\varphi }}}}_{\mathrm{1}}^{*}+{\displaystyle \frac{{\partial }^{\mathrm{2}}{\stackrel{\mathrm{‾}}{\stackrel{\mathrm{̃}}{\mathit{\varphi }}}}_{\mathrm{1}}^{*}}{\partial y^{*\mathrm{2}}}}={\displaystyle }& {\displaystyle \frac{\mathrm{1}}{p}}{\mathit{\theta }}_{\mathrm{1}}^{\mathrm{2}}{\mathit{\lambda }}_i(-\mathrm{1}{)}^i\left(h_{\mathrm{31}}^{*}+h_{\mathrm{23}}^{*}{y}^{*}\right)-\sum _{m=\mathrm{1}}^{\mathrm{\infty }}{\mathrm{\Delta }}_{\mathrm{1}}\left[C_{\mathrm{1}m}{E}_{\mathrm{1}}\left(m,y^{*}\right)\right.\\ {\displaystyle }& {\displaystyle }\left.+F_{\mathrm{1}}\left(m,y^{*}\right)\right]\left({\displaystyle \frac{\mathrm{1}}{\mathrm{2}}}-{\displaystyle \frac{\sin \mathrm{2}{\mathit{\lambda }}_i}{\mathrm{4}{\mathit{\lambda }}_i}}\right)\\ {\displaystyle }& {\displaystyle }-{\displaystyle \frac{\mathrm{1}}{p}}\sum _{k=\mathrm{1}}^M{Q}_{\mathrm{1}k}^{*}\sin \left({\mathit{\lambda }}_i{x}_{\mathrm{1}k}^{*}\right)\mathit{\delta }\left(y^{*}-y_{\mathrm{1}k}^{*}\right),\end{array}$$

where ${\mathit{\mu }}_i=\sqrt{{\mathit{\theta }}_{\mathrm{1}}^{\mathrm{2}}{\mathit{\lambda }}_i^{\mathrm{2}}+p}$, with the Laplace transform of ${\stackrel{\mathrm{‾}}{\mathit{\varphi }}}_{\mathrm{1}}^{*}$ defined as

$$\begin{array}{}\text{(B7)}& {\displaystyle }{\stackrel{\mathrm{‾}}{\stackrel{\mathrm{̃}}{\mathit{\varphi }}}}_{\mathrm{1}}^{*}\left(i,y^{*},p\right)=\underset{\mathrm{0}}{\overset{\mathrm{\infty }}{\int }}{\stackrel{\mathrm{‾}}{\mathit{\varphi }}}_{\mathrm{1}}^{*}\left(i,y^{*},t\right)e^{-p{t}^{*}}\text{d}{t}^{*}.\end{array}$$

Similarly, the Laplace transform of Eq. (B3) is obtained as

$$\begin{array}{ll}\text{(B8)}& {\displaystyle }& {\displaystyle }-{\mathrm{\Omega }}_{\mathrm{2}j}^{\mathrm{2}}{\stackrel{\mathrm{‾}}{\mathit{\varphi }}}_{\mathrm{2}}^{*}+{\displaystyle \frac{\mathrm{1}}{p}}{\mathit{\theta }}_{\mathrm{2}}^{\mathrm{2}}{\mathit{\alpha }}_j{\left(-\mathrm{1}\right)}^j\left(H_{\mathrm{31}}^{*}+H_{\mathrm{23}}^{*}{y}^{*}\right)+{\displaystyle \frac{{\partial }^{\mathrm{2}}{\stackrel{\mathrm{‾}}{\mathit{\varphi }}}_{\mathrm{2}}^{*}}{\partial y^{*\mathrm{2}}}}={\displaystyle }& {\displaystyle \frac{K_{y\mathrm{1}}{S}_{s\mathrm{2}}}{K_{y\mathrm{2}}{S}_{s\mathrm{1}}}}(p{\stackrel{\mathrm{‾}}{\stackrel{\mathrm{̃}}{\mathit{\varphi }}}}_{\mathrm{2}}^{*}-{\stackrel{\mathrm{‾}}{\mathit{\varphi }}}_{\mathrm{2}s}^{*})-{\displaystyle \frac{\mathrm{1}}{p}}\sum _{l=\mathrm{1}}^N{Q}_{\mathrm{2}l}^{*}\cos \left({\mathit{\alpha }}_j{x}_{\mathrm{2}l}^{*}\right)\mathit{\delta }\left(y^{*}-y_{\mathrm{2}l}^{*}\right),\end{array}$$

where ${\stackrel{\mathrm{‾}}{\mathit{\varphi }}}_{\mathrm{2}s}^{*}$ is the steady-state solution of region 2. Thus, Eq. (B8) can be written as

$$\begin{array}{ll}\text{(B9)}& {\displaystyle }& {\displaystyle }-{\mathit{\theta }}_j^{\mathrm{2}}{\stackrel{\mathrm{‾}}{\stackrel{\mathrm{̃}}{\mathit{\varphi }}}}_{\mathrm{2}}^{*}+{\displaystyle \frac{{\partial }^{\mathrm{2}}{\stackrel{\mathrm{‾}}{\stackrel{\mathrm{̃}}{\mathit{\varphi }}}}_{\mathrm{2}}^{*}}{\partial y^{*\mathrm{2}}}}={\displaystyle }& {\displaystyle }-{\displaystyle \frac{\mathrm{1}}{p}}{\mathit{\theta }}_{\mathrm{2}}^{\mathrm{2}}{\mathit{\alpha }}_j(-\mathrm{1}{)}^j\left(H_{\mathrm{31}}^{*}+H_{\mathrm{23}}^{*}{y}^{*}\right)-{\displaystyle \frac{K_{y\mathrm{1}}{S}_{s\mathrm{2}}}{K_{y\mathrm{2}}{S}_{s\mathrm{1}}}}\sum _{n=\mathrm{1}}^{\mathrm{\infty }}{\mathrm{\Delta }}_{\mathrm{2}}\\ {\displaystyle }& {\displaystyle }\left[D_{\mathrm{2}n}{E}_{\mathrm{2}}\left(m,y^{*}\right)+F_{\mathrm{2}}\left(m,y^{*}\right)\right]\cos \left[{\mathit{\alpha }}_j\left(x^{*}-l_{\mathrm{2}}^{*}\right)\right]\\ {\displaystyle }& {\displaystyle }-{\displaystyle \frac{\mathrm{1}}{p}}\sum _{l=\mathrm{1}}^N{Q}_{\mathrm{2}l}^{*}\cos \left({\mathit{\alpha }}_j{x}_{\mathrm{2}l}^{*}\right)\mathit{\delta }\left(y^{*}-y_{\mathrm{2}l}^{*}\right),\end{array}$$

where ${\mathit{\theta }}_j=\sqrt{{\mathit{\theta }}_{\mathrm{2}}^{\mathrm{2}}{\mathit{\alpha }}_j^{\mathrm{2}}+p{s}_{s\mathrm{2}}{k}_{y\mathrm{1}}/s_{s\mathrm{1}}{k}_{y\mathrm{2}}}$, with the Laplace transform of ${\stackrel{\mathrm{‾}}{\mathit{\varphi }}}_{\mathrm{2}}^{*}$ defined as

$$\begin{array}{}\text{(B10)}& {\displaystyle }{\stackrel{\mathrm{‾}}{\stackrel{\mathrm{̃}}{\mathit{\varphi }}}}_{\mathrm{2}}^{*}\left(j,y^{*},p\right)=\underset{\mathrm{0}}{\overset{\mathrm{\infty }}{\int }}{\stackrel{\mathrm{‾}}{\mathit{\varphi }}}_{\mathrm{2}}^{*}\left(j,y^{*},t\right)e^{-p{t}^{*}}\text{d}{t}^{*}.\end{array}$$

The general solution of Eq. (B6) can be expressed as

$$\begin{array}{}\text{(B11)}& {\displaystyle }{\stackrel{\mathrm{‾}}{\stackrel{\mathrm{̃}}{\mathit{\varphi }}}}_{\mathrm{1}}^{*}\left(i,y^{*},p\right)=T_{\mathrm{1}i}{e}^{{\mathit{\mu }}_i{y}^{*}}+T_{\mathrm{2}i}{e}^{-{\mathit{\mu }}_i{y}^{*}}+{\stackrel{\mathrm{‾}}{\stackrel{\mathrm{̃}}{\mathit{\varphi }}}}_{\mathrm{1}p}^{*}\left(i,y^{*},p\right),\end{array}$$

where the particular solution ${\stackrel{\mathrm{‾}}{\stackrel{\mathrm{̃}}{\mathit{\varphi }}}}_{\mathrm{1}p}^{*}\left(i,y^{*},p\right)$ is

$$\begin{array}{ll}\text{(B12)}& {\displaystyle }& {\displaystyle }{\stackrel{\mathrm{‾}}{\stackrel{\mathrm{̃}}{\mathit{\varphi }}}}_{\mathrm{1}p}^{*}\left(i,y^{*},p\right)={\displaystyle }& {\displaystyle \frac{e^{{\mathit{\mu }}_i{y}^{*}}}{\mathrm{2}{\mathit{\mu }}_i}}\int e^{-{\mathit{\mu }}_i{y}^{*}}{\mathrm{\Delta }}_{\mathrm{1}y}\left(i,y^{*},p\right)\text{d}{y}^{*}\\ {\displaystyle }& {\displaystyle }-{\displaystyle \frac{e^{-{\mathit{\mu }}_i{y}^{*}}}{\mathrm{2}{\mathit{\mu }}_i}}\int e^{{\mathit{\mu }}_i{y}^{*}}{\mathrm{\Delta }}_{\mathrm{1}y}\left(i,y^{*},p\right)\text{d}{y}^{*},\end{array}$$

with

$$\begin{array}{ll}\text{(B13)}& {\displaystyle }& {\displaystyle }{\mathrm{\Delta }}_{\mathrm{1}y}\left(i,y^{*},p\right)={\displaystyle }& {\displaystyle \frac{\mathrm{1}}{p}}{\mathit{\theta }}_{\mathrm{1}}^{\mathrm{2}}{\mathit{\lambda }}_i(-\mathrm{1}{)}^i\left(h_{\mathrm{31}}^{*}+h_{\mathrm{23}}^{*}{y}^{*}\right)-\sum _{m=\mathrm{1}}^{\mathrm{\infty }}{\mathrm{\Delta }}_{\mathrm{1}}\left[C_{\mathrm{1}m}{E}_{\mathrm{1}}\left(m,y^{*}\right)\right.\\ {\displaystyle }& {\displaystyle }\left.+F_{\mathrm{1}}\left(m,y^{*}\right)\right]\left({\displaystyle \frac{\mathrm{1}}{\mathrm{2}}}-{\displaystyle \frac{\sin \mathrm{2}{\mathit{\lambda }}_i}{\mathrm{4}{\mathit{\lambda }}_i}}\right)\\ {\displaystyle }& {\displaystyle }-{\displaystyle \frac{\mathrm{1}}{p}}\sum _{k=\mathrm{1}}^M{Q}_{\mathrm{1}k}^{*}\sin \left({\mathit{\lambda }}_i{x}_{\mathrm{1}k}^{*}\right)\mathit{\delta }\left(y^{*}-y_{\mathrm{1}k}^{*}\right).\end{array}$$

Moreover, Eq. (B9) can also be expressed as

$$\begin{array}{}\text{(B14)}& {\displaystyle }{\stackrel{\mathrm{‾}}{\stackrel{\mathrm{̃}}{\mathit{\varphi }}}}_{\mathrm{2}}^{*}\left(j,y^{*},p\right)=T_{\mathrm{1}j}{e}^{{\mathit{\theta }}_j{y}^{*}}+T_{\mathrm{2}j}{e}^{-{\mathit{\theta }}_j{y}^{*}}+{\stackrel{\mathrm{‾}}{\stackrel{\mathrm{̃}}{\mathit{\varphi }}}}_{\mathrm{2}p}^{*}\left(j,y^{*},p\right),\end{array}$$

in which ${\stackrel{\mathrm{‾}}{\stackrel{\mathrm{̃}}{\mathit{\varphi }}}}_{\mathrm{2}p}^{*}\left(j,y^{*},p\right)$ is

$$\begin{array}{ll}\text{(B15)}& {\displaystyle }& {\displaystyle }{\stackrel{\mathrm{‾}}{\stackrel{\mathrm{̃}}{\mathit{\varphi }}}}_{\mathrm{2}p}^{*}\left(j,y^{*},p\right)={\displaystyle }& {\displaystyle \frac{e^{{\mathit{\theta }}_j{y}^{*}}}{\mathrm{2}{\mathit{\theta }}_j}}\int e^{-{\mathit{\theta }}_j{y}^{*}}{\mathrm{\Delta }}_{\mathrm{2}y}\left(j,y^{*},p\right)\text{d}{y}^{*}\\ {\displaystyle }& {\displaystyle }-{\displaystyle \frac{e^{-{\mathit{\theta }}_j{y}^{*}}}{\mathrm{2}{\mathit{\theta }}_j}}\int e^{{\mathit{\theta }}_j{y}^{*}}{\mathrm{\Delta }}_{\mathrm{2}y}\left(i,y^{*},p\right)\text{d}{y}^{*},\end{array}$$

with

$$\begin{array}{ll}\text{(B16)}& {\displaystyle }& {\displaystyle }{\mathrm{\Delta }}_{\mathrm{2}y}\left(j,y^{*},p\right)={\displaystyle }& {\displaystyle }-{\displaystyle \frac{\mathrm{1}}{p}}{\mathit{\theta }}_{\mathrm{2}}^{\mathrm{2}}{\mathit{\alpha }}_j{\left(-\mathrm{1}\right)}^j\left(H_{\mathrm{31}}^{*}+H_{\mathrm{23}}^{*}{y}^{*}\right)\\ {\displaystyle }& {\displaystyle }-{\displaystyle \frac{K_{y\mathrm{1}}{S}_{s\mathrm{2}}}{K_{y\mathrm{2}}{S}_{s\mathrm{1}}}}\sum _{n=\mathrm{1}}^{\mathrm{\infty }}{\mathrm{\Delta }}_{\mathrm{2}}\left[D_{\mathrm{2}n}{E}_{\mathrm{2}}\left(m,y^{*}\right)+F_{\mathrm{2}}\left(m,y^{*}\right)\right]\\ {\displaystyle }& {\displaystyle }\cos \left[{\mathit{\alpha }}_j\left(x^{*}-l_{\mathrm{2}}^{*}\right)\right]-{\displaystyle \frac{\mathrm{1}}{p}}\sum _{l=\mathrm{1}}^N{Q}_{\mathrm{2}l}^{*}\cos \left({\mathit{\alpha }}_j{x}_{\mathrm{2}l}^{*}\right)\mathit{\delta }\left(y^{*}-y_{\mathrm{2}l}^{*}\right).\end{array}$$

On the basis of Eq. (A5), the coefficient T_2i in Eq. (B11) can be determined in terms of T_1i as

$$\begin{array}{ll}\text{(B17)}& {\displaystyle }& {\displaystyle }{T}_{\mathrm{2}i}={\displaystyle }& {\displaystyle }-h_{\mathrm{31}}^{*}{\displaystyle \frac{{\left(-\mathrm{1}\right)}^i}{{\mathit{\lambda }}_i}}+{\displaystyle \frac{h_{\mathrm{31}}^{*}}{p{\mathit{\mu }}_i^{\mathrm{2}}}}\left[{\mathit{\theta }}_{\mathrm{1}}^{\mathrm{2}}{\mathit{\lambda }}_i{\left(-\mathrm{1}\right)}^i\right]\\ {\displaystyle }& {\displaystyle }+{\mathrm{\Delta }}_{\mathrm{1}}{\left(-\mathrm{1}\right)}^i{C}_{\mathrm{1}m}\left({\displaystyle \frac{\mathrm{1}}{\mathrm{2}}}-{\displaystyle \frac{\sin \mathrm{2}{\mathit{\lambda }}_i}{\mathrm{4}{\mathit{\lambda }}_i}}\right)-T_{\mathrm{1}i}.\end{array}$$

The solution for hydraulic head distribution in region 1 is given as Eq. (26), which is obtained by substituting Eq. (B17) into Eq. (B11) and then taking the following inverse Fourier transform to Eq. (B11) denoted as

$$\begin{array}{}\text{(B18)}& {\displaystyle }{\stackrel{\mathrm{̃}}{\mathit{\varphi }}}_{\mathrm{1}}^{*}\left(x^{*},y^{*},p\right)={\sum }_{i=\mathrm{0}}^{\mathrm{\infty }}{\stackrel{\mathrm{‾}}{\stackrel{\mathrm{̃}}{\mathit{\varphi }}}}_{\mathrm{1}}^{*}\left(i,y^{*},p\right)\sin \left({\mathit{\lambda }}_i{x}^{*}\right),\end{array}$$

with

$$\begin{array}{ll}\text{(B19)}& {\displaystyle }& {\displaystyle }{w}_{\mathrm{1}i}^{*}={\displaystyle }& {\displaystyle }{\left.T_{\mathrm{1}i}(e^{{\mathit{\mu }}_i{y}^{*}}-\mathrm{1})\right|}_{y^{*}=d_{\mathrm{1}}^{*}}-{\displaystyle \frac{\mathrm{1}}{\mathrm{2}{\mathit{\mu }}_{i}p}}\sum _{k=\mathrm{1}}^M{Q}_{\mathrm{1}k}^{*}\sin \left({\mathit{\lambda }}_i{x}_{\mathrm{1}k}^{*}\right)\\ {\displaystyle }& {\displaystyle }{\left.\left[e^{{\mathit{\mu }}_i(y^{*}-y_{\mathrm{1}k}^{*})}-e^{{\mathit{\mu }}_i(y_{\mathrm{1}k}^{*}-y^{*})}\right]\right|}_{y^{*}=d_{\mathrm{1}}^{*}}.\end{array}$$

Similarly, T_1j in Eq. (B14) can be obtained based on Eq. (A9) as

$$\begin{array}{ll}\text{(B20)}& {\displaystyle }& {\displaystyle }{T}_{\mathrm{1}j}={\displaystyle }& {\displaystyle }{T}_{\mathrm{2}j}{e}^{-\mathrm{2}{\mathit{\theta }}_j}-\left[{\displaystyle \frac{{\mathit{\theta }}_{\mathrm{2}}^{\mathrm{2}}{\mathit{\alpha }}_j{\left(-\mathrm{1}\right)}^j}{p}}+{\displaystyle \frac{K_{y\mathrm{1}}{S}_{s\mathrm{2}}{\left(-\mathrm{1}\right)}^j}{K_{y\mathrm{2}}{S}_{s\mathrm{1}}{\mathit{\alpha }}_j}}\right]{\displaystyle \frac{H_{\mathrm{21}}^{*}{e}^{-{\mathit{\theta }}_j}}{{\mathit{\theta }}_j^{\mathrm{2}}}}\\ {\displaystyle }& {\displaystyle }+{\displaystyle \frac{H_{\mathrm{21}}^{*}{e}^{-{\mathit{\theta }}_j}{\left(-\mathrm{1}\right)}^j}{{\mathit{\alpha }}_{j}p}}+{\displaystyle \frac{\mathrm{1}}{p}}\sum _{l=\mathrm{1}}^N{Q}_{\mathrm{2}l}^{*}\cos \left({\mathit{\alpha }}_j{x}_{\mathrm{2}l}^{*}\right){\displaystyle \frac{e^{-{\mathit{\theta }}_j{y}_{\mathrm{2}l}^{*}}-e^{{\mathit{\theta }}_j(y_{\mathrm{2}l}^{*}-\mathrm{2})}}{\mathrm{2}{\mathit{\theta }}_j}}.\end{array}$$

The solution for region 2 is Eq. (27), which is acquired by substituting Eq. (B20) into Eq. (B14) and then taking the following inverse Fourier transform to Eq. (B14) expressed as

$$\begin{array}{}\text{(B21)}& {\displaystyle }{\stackrel{\mathrm{̃}}{\mathit{\varphi }}}_{\mathrm{2}}^{*}\left(x^{*},y^{*},p\right)={\sum }_{j=\mathrm{0}}^{\mathrm{\infty }}{\stackrel{\mathrm{‾}}{\stackrel{\mathrm{̃}}{\mathit{\varphi }}}}_{\mathrm{2}}^{*}\left(j,y^{*},p\right)\cos \left({\mathit{\alpha }}_j{x}^{*}\right),\end{array}$$

with

$$\begin{array}{}\text{(B22)}& {\displaystyle }{w}_{\mathrm{2}j}^{*}=T_{\mathrm{2}j}.\end{array}$$

Furthermore, the coefficients of $w_{\mathrm{1}i}^{*}$ and $w_{\mathrm{2}j}^{*}$ can be simultaneously determined by Eqs. (A10) and (A11). The results are respectively given in Eqs. (36) and (37).