Analytical Solutions for Steady-State Multiwell Aquifer Tests in Rectangular Aquifers by Using Double Fourier Transform: A Case Study in the Ordos Plateau, China

Straightforward solutions have long been expected for the analysis of multiwell aquifer tests. In this paper, we derive series analytical solutions of steady-state groundwater flow in a rectangular-shaped aquifer with pumping/injection wells for both confined and unconfined conditions. Double Fourier Transform (DFT) technique is applied to deal with different combinations of impermeable and specified head boundaries on sides. The obtained solutions are compact and concise in mathematics and flexible in terms of well number, well locations, and pumping/injection rates. Hatoucaidang, a groundwater resource field in the Ordos Plateau, Northwestern China, is introduced as a field case study, where a multiwell aquifer test was conducted. One of the analytical solutions derived herein is used to estimate hydraulic conductivities by applying a direct calculation method and a least square estimation method regarding observed versus calculated drawdowns. By comparing with nearby single-well pumping tests, the reliability of the derived analytical solutions is proven. This study facilitates utilizing the multiwell aquifer test to analyze the general behavior of groundwater movement in aquifer systems.


Introduction
Aquifer tests are widely adopted as a standard practice in groundwater science and engineering, which can be used not only to characterize the general behavior of groundwater movement in aquifer systems but also to more precisely estimate hydrogeological parameters of the study site [1][2][3][4]. In many cases, a single-well aquifer test (i.e., one pumping/injection well with/without observation wells) can basically meet the demands of the conventional hydrogeological surveys or research. However, with the development of society, requirements of the accuracy in hydrogeological survey or research are getting higher. As one of the most efficient and promising methods, multiwell aquifer tests are increasingly needed. A multiwell aquifer test is usually characterized by a high rate of pumping/injection, a deep drawdown/uplift, and long-term operations. Thus, it has obvious advantages in solving such problems as identifying the source of groundwater recharge and the principle direction of regional groundwater flow, estimating hydrogeological parameters, analyzing hydraulic connections between groundwater and surface water, designing the framework for pump-and-treat, etc.(e.g., [5][6][7][8][9]).
Compared with a single-well test, multiwell tests are more likely to be influenced by the bounded hydrogeological boundaries such as impervious rocks, fault zones, perennial rivers, and coastlines [10,11] and by the mutual interference of multiple pumping/injection wells. By using the image well approach and conformal mapping approach, theoretically, analytical solutions of a single well can be transformed to fit for the multiwell aquifer tests (e.g., [10,[12][13][14][15]). However, as well number increases and pumping/injection wells are arbitrarily located, in theory, an infinite number of image wells would be required to ensure that the solution accuracy and the application of corresponding analytical solutions would become toughly complex [10,16]. With the recent advance of computer technology, numerical modeling serves the survey and research of hydrogeology intensely, and the hydrogeological parameters of aquifers can be obtained by inversion of a groundwater numerical model. The numerical model is quite versatile for parameter estimation using a calibration process, but it also has the disadvantages of unknown numerical errors, nonuniqueness of identified parameter values, and sometimes convergence issues. Therefore, it is of theoretical and practical value to derive straightforward analytical solutions with compactness, conciseness, and clear physical meanings for the problem.
As a standard and well-known practice, Double Fourier Transform (DFT) has been widely used in the field of hydrogeology [8,17,18]. In this study, we apply the DFT method to derive series analytical solutions of steady-state groundwater flow in rectangular aquifers in the presence of multiple wells. The series analytical solutions can, respectively, deal with five hydrogeological scenarios with different combinations of impermeable and specified head boundaries and can be conveniently used to estimate hydraulic conductivities, which has been applied in the designation of a drinking groundwater resource field in the Ordos Plateau, Northwestern China.

Methods
The issue of pumping in rectangular aquifers has been widely studied both in the confined and unconfined aquifers previously (e.g., [12,16,[19][20][21]). However, these solutions are mainly developed for a single well; effects of multiwell pumping/injections on the groundwater flow were not considered, especially for cases that the pumping/injection wells are arbitrarily located and at different rates.
Here, we consider that multiple fully penetrated pumping/injection wells are located within a confined (or unconfined) aquifer of rectangular shape in horizontal, as shown in Figure 1(a). The x and y axes are aligned along the lower and left boundaries, respectively, with the origin at the lower left corner. The length (in the x-direction) and width (in the y-direction) of the aquifer are L x and L y , respectively. A pumping/injection well is settled at an arbitrary location with coordinates of ðx i , y i Þ and operates at any reasonable constant rate of Q i . Groundwater flow is assumed at the steady-state situation. For simplification, surface recharge/discharge and regional groundwater flow are neglected in this study.
Consistent with former studies [12,16,[19][20][21], five hydrogeological scenarios with different combinations of impermeable and specified head boundaries are discussed. Specifically, in scenario 1 (Figure 1(b)), the aquifer is bounded by four specified head boundaries; in scenario 2 ( Figure 1(c)), the aquifer is bounded by three specified head boundaries and one impermeable boundary; in scenario 3 ( Figure 1(d)), the aquifer is bounded by two parallel impermeable boundaries and two parallel specified head boundaries; in scenario 4 ( Figure 1(e)), the aquifer is bounded by two pairs of orthogonal impermeable and specified head boundaries; in scenario 5 (Figure 1(f)), the aquifer is bounded by one specified head boundary and three impermeable boundaries. The scenario that an aquifer fully bounded by four impermeable boundaries is not considered, because it will never lead to a steady-state condition [12].

Solutions for a Confined
Aquifer. The governing equation of drawdown induced by multiple wells in a confined aquifer can be given as where T x and T y are the hydraulic transmissivities in the x-and y-directions, respectively, s is the drawdown, w gives the number of wells, Q i refers to the pumping rate of the ith well, δ is the Dirac delta function, and x i and y i stand for the coordinate of the ith well. Note that Q i < 0 stands for injection and Q i > 0 for pumping. As the governing equation is a second-order linear partial differential equation, the superposition principle can be used. For each hydrogeological scenario, by using methods of Separation of Variables (SOF) and DFT, the analytical solution could be obtained. A detailed derivation process can be found in the appendix.
In order to reveal intrinsic factors controlling the general behaviors of groundwater flow, we introduce the following dimensionless parameters: and the series multiwell analytical solutions could be expressed as where Wðx i , y i , x, y, α, βÞ is the well function of x i , y i , x, y, α, and β. Specifically, in scenario 1, the well function can be expanded as In scenario 2, the well function becomes In scenario 3, the well function is given as where δ m = 2 for m = 0 and δ m = 1 for m ≠ 0, respectively. Note that m starts from zero in this situation. In scenario 4, the well function becomes In scenario 5, the well function is given as where δ n = 2 for n = 0 and δ n = 1 for n ≠ 0, respectively. Note that n starts from zero in this situation.
In the above solutions, the mathematical appearances are generally the same. The difference mainly exists in the well functions. Well functions essentially are the expression of characteristic functions of the boundary conditions. Therefore, these series analytical solutions have consistency in the form of mathematics, and one can be easily modified to another by using the corresponding characteristic functions according to the boundary conditions.

Geofluids
It should be mentioned that Equations (3) and (4) for scenario 1 are identical to those given by Yeo and Lee [8] since they deal with the same boundary conditions. In scenario 5, although it is a cross-sectional model in Wang et al. [17] while a planar model in this paper, the boundary conditions of the two are the same. Thus, the solution of using Equations (3) and (8) is essentially identical to that given by Wang et al. [17]. For the conditions with impermeable boundaries in scenarios 3 and 4, Yeo and Lee [8] recommended using the image well approach, but we can give the straightforward solutions with Equations (6) and (7).

Solutions for an Unconfined
Aquifer. The governing equation of water table height in an unconfined aquifer in the presence of multiple wells can be given as where K x and K y are the hydraulic conductivities in the xand y-directions and h is the water table height above the flat aquifer bottom.

Geofluids
Equation (9) is a second-order nonlinear partial differential equation but can be linearized easily by introducing f = h 2 /2 as follows: which is similar to that shown in Equation (1). Therefore, by properly modifying the multiwell analytical solutions of a confined aquifer, corresponding analytical solutions of an unconfined aquifer can be obtained: where h 0 is the initial water table height (a constant), and the well function Wðx i , y i , x, y, α, βÞ can be determined with Equations (3)-(8) for different scenarios as that defined for the study of a confined aquifer. Note that Equation (11) is not a straightforward solution for the drawdown s. But in a thick unconfined aquifer, it can be generally simplified, when the change in the water table is significantly small in comparison with the initial height, i.e., s ≪ h 0 , we could use Δf ≈ h 0 s, which leads to It should be noted that in this study, the transient behaviors of groundwater flow triggered by aquifer tests are not considered. The steady-state solutions are fine, but it means that valuable transient data sets collected during aquifer tests cannot be used. Therefore, no data of storativity, specific yield, or aquifer diffusivity can be obtained. These parameters are vital for the management of groundwater resources, the evaluation of surface water-groundwater interactions, and the protection of the ecological environment [22][23][24][25]. Solutions regarding transient behaviors should be considered in the future.

A Field Case Study
3.1. Study Area. The Ordos Plateau is located in Northwestern China (Figure 2(a)). Due to the arid to semiarid climate, groundwater is the only useable water resource in most parts [26]. In recent years, to satisfy the increasing demand for drinking water, the local government began a project of  5 Geofluids pumping groundwater from a site named Hatoucaidang with a group of wells (Figure 2(b)).
The groundwater resource field in Hatoucaidang, with a total area of about 1300 km 2 , is located in the hinterland of the Mu Us Desert of the Ordos Plateau, southwest of the Ordos city (Figure 2(c)). The landscape is characterized by sparse vegetation and relatively flat topography. The target exploitation aquifer is formed by the Quaternary sediments where groundwater depth is generally less than 3 m but the depth of the aquifer bottom is larger than 100 m. These Quaternary sediments act as a porous unconfined aquifer with lithology dominated by medium to fine sands. Underlying this unconfined aquifer, the Cretaceous sandstones exist, which have relatively very low permeability in comparison with the Quaternary sediments. Thus, the top of sandstones could be plausibly regarded as the impervious boundary of the unconfined aquifer.

Multiwell Aquifer Test.
In order to ascertain the potential yield capacity of the groundwater resource field, a multiwell aquifer test was scheduled. The test site was located on the east side of the proposed groundwater resource field (Figure 2(c)). There were 5 pumping wells (labeled as HT07, HT08, HT09, HT10, and HD27) and 18 observation wells included in the test. Locations of each well are shown in Figure 3. The observation wells can be categorized into two types: piezometric wells (labeled with Obs in front) drilled in this test and observed with the automatic water level gauge and local domestic wells (labeled with MJ in front) formerly drilled by the residents and observed manually. During the pumping process, flow rates of each pumping well were kept constant. The multiwell aquifer test lasts for about 20 days with a total abstraction of about 360000 m 3 .

Application of the Analytical Solution.
The hydraulic conductivity of the aquifer could be estimated by the analytical solution derived in the former section. Assuming that the Quaternary sediments are homogeneous and isotropic, we have α = 1 and K = K x = K y . The average aquifer thickness h 0 is 124.6 m according to the drilling logs of each pumping well ( Table 1). The rectangular calculating area is centered at the pumping wells and extends to a certain distance in the x-or y-direction. The extending distance is determined by checking the drawdowns in the domestic observation wells far away from the pumping site ( Figure 3, Table 2). Drawdown less than 0.1 m is considered to be not induced by pumping but caused by some natural processes such as seasonal dynamics. Therefore, we assume a model with the length L x = 4500 m and the width L y = 3000 m, i.e., β = 2/3 ( Figure 3).
Here, we choose to use the analytical solution regarding four specified head boundaries in an unconfined aquifer (Equations (11) and (4)). In this case, the Quaternary unconfined aquifer is thick, and the change in water table height is relatively small compared to the thickness of the aquifer. Accordingly, we can use the simplified formula, Equation (12), to directly calculate the hydraulic conductivity as follows: where the drawdown s of a selected observation well   (Table 1), the coordinates of each pumping well ðx i , y i Þ and a selected observation well ð x, yÞ, and the drawdown of a selected observation well s (Figure 3, Table 2) into Equation (13), the hydraulic conductivity can be estimated (the domestic wells were not used because of uncertainty in observations). As shown in Table 2, this results in 9 different hydraulic conductivities of K values in the range of 25.84~59.27 m/d due to the nature of heterogeneity and noise in the observations. To average, the K value is 38.23 m/d. Another way is introducing a single effective value of the hydraulic conductivity in Equation (12) to calculate the drawdown values for all the wells that reach a best fitness with the observations. Here, the sum of squares of errors (SSE) between the observed and calculated drawdowns is chosen as the objective function so that the method of least squares can be applied to determine the optimized value of the hydraulic conductivity: where n is the number of observation wells, obs is the abbreviation of observation, and cal is the abbreviation of calculation. As a result, the best fitting hydraulic conductivity of the K value is 31.51 m/d. The spatial distribution of drawdown calculated by the multiwell analytical solution with the best fitting K value of 31.51 m/d is shown in Figure 3. Due to the nature of heterogeneity, the observed drawdowns are not exactly the same as the calculated ones, but the two groups have similar values with a maximum difference of 0.33 m (Figure 4). Note that the best fitting K value given by the least square method is not significantly different from the average of K values by a direct calculation method with respect to Equation (13) (38.23 m/d) but may be more reasonable because it matches the observed drawdown in an integrated way.

Discussion
Before the multiwell aquifer test, 4 sets of single-well aquifer tests were conducted nearby. Each single-well aquifer test included one pumping well and three observation wells. Key records of each set of the single-well aquifer test are shown in (Table 3). The hydraulic conductivity was inversely estimated with the Thiem equation of an unconfined aquifer [3,4] from the final steady-state drawdown, leading to K values in the range of 13.13~27.71 m/d with the average point at 19.17 m/d. Note that the hydraulic conductivity seems to be underestimated by these singlewell aquifer tests in comparison with the results obtained from the multiwell aquifer test. This underestimation may be caused by the significantly limited influence zone of a single pumping well. The multiwell aquifer test provided a more reasonable estimation of the hydraulic conductivity in a larger influenced area.

Conclusions
Analytical solutions for multiwell aquifer tests are of important theoretical significance and practical application values. In this study, by using the Double Fourier Transform (DFT), we derived series analytical solutions of steady-state groundwater flow in a rectangular aquifer in the presence of multiple pumping or injection wells. These series multiwell analytical solutions can deal with five hydrological scenarios of different combinations of impermeable and specified head boundaries in both the confined and unconfined aquifers and have no restrictions on well number, well locations, and pumping/injection rates. Moreover, they are relatively compact and concise in mathematics compared to previous studies using the image well approach. Hatoucaidang, a groundwater resource field in the Ordos Plateau, Northwestern China, is introduced as a field case study. A multiwell aquifer test with 5 pumping wells and 9 observation wells was conducted. The analytical solution regarding four specified head boundaries in an unconfined aquifer (i.e., Equations (13)) is used to directly calculate hydraulic conductivities for each selected observation well. In addition, the least square method is introduced to estimate a single effective value of the hydraulic conductivity using the drawdown values for all the wells that reach the best fitness with the observations. As a result, the hydraulic conductivity by the least square method is similar to those by the direct calculation method. Moreover, 4 sets of single-well aquifer tests conducted nearby were used to prove the reliability of the derived analytical 7 Geofluids solutions. The hydraulic conductivities given by these single-well aquifer tests are similar but smaller in comparison with the results obtained from the multiwell aquifer test, which may be due to the significantly limited influence zone of a single pumping well.
This study facilitates utilizing the multiwell aquifer test to analyze the general behavior of groundwater movement in aquifer systems.

Appendix A Detailed Derivation of Multiwell Analytical Solutions
In scenario 1 (Figure 1(b)), the mathematical model of the multiwell aquifer test is ðA:10Þ Thus, the governing equation can be transformed into ðA:11Þ After rearrangement, we can obtain ðA:14Þ By putting the following dimensionless parameters into 8 Geofluids Equation (A.14), the dimensionless solution of drawdown s is where subscript C denotes confined aquifer and S1 denotes scenario 1.
In scenario 2, the aquifer is bounded by one impermeable and three specified head boundaries (Figure 1(c)). According to the boundary conditions, for variable x, the corresponding characteristic function is ðA:18aÞ and for variable y, the corresponding characteristic function is ðA:18bÞ Using the same techniques as scenario 1, the analytical solution of drawdown s could be derived as Þπ y 2 ,
In scenario 3, the aquifer is bounded by two parallel impermeable boundaries in the x-direction and two parallel specified head boundaries in the y-direction (Figure 1(d)). According to the boundary conditions, for variable x, the corresponding characteristic function is ðA:21aÞ and for variable y, the corresponding characteristic function is ðA:21bÞ Using the same techniques as scenario 1, the analytical solution of drawdown s could be derived as

ðA:23Þ
where subscript C denotes confined aquifer and S3 denotes scenario 3. Note that m starts from zero in W CS3 ðx i , y i , x, y, α, βÞ.
In scenario 4, the aquifer is bounded by two pairs of orthotropic impermeable and orthotropic specified head boundaries (Figure 1(e)). According to the boundary conditions, for variable x, the corresponding characteristic function is ðA:24bÞ Using the same techniques as scenario 1, the analytical solution of drawdown s could be derived as where subscript C denotes confined aquifer and S4 denotes scenario 4. In scenario 5, the aquifer is bounded by one specified head boundary and three impermeable boundaries (Figure 1(f)). According to the boundary conditions, for variable x, the corresponding characteristic function is ðA:27bÞ Using the same techniques as scenario 1, the analytical solution of drawdown s could be derived as where subscript C denotes confined aquifer and S5 denotes scenario 5. Note that n starts from zero in W CS5 ðx i , y i , x, y, α, βÞ.
In the above solutions, the mathematical appearances are generally the same. The difference mainly exists in the well functions. Well functions essentially are the expression of characteristic functions of the boundary conditions. Therefore, these series analytical solutions can be expressed in the general form where Wðx i , y i , x, y, α, βÞ is the well function as a function of different boundary combinations.

Data Availability
The field data regarding the multiwell aquifer test in the Ordos Plateau are generally included within the paper. More detailed information is available from the authors upon request.

Conflicts of Interest
The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.