An Analytical Solution for Non-Darcian Flow on a Constant Head Packer Test in the Interlayer Staggered Zone

,


Introduction
The permeability of rock mass is mainly determined by the development of fractures, and the packer test has been widely used to characterize the permeability of a lowpermeability aquifer. The slug test and constant head injection test [1] are two common in situ packer tests to estimate the permeability of fractured rock conveniently. The constant head packer test is commonly applied for lowpermeability aquifers, and the non-Darcian flow is prone to occur with high gradient and high velocity during the test.
In general, the linear relationship between the velocity of groundwater flow and the hydraulic gradient is called Darcy's law. When the velocity is too high or too low, however, the non-Darcian flow may occur. According to the flow velocity, there are two types of non-Darcian flow: the prelinear flow and the postlinear flow. The prelinear flow occurs when the flow at low velocity, especially in clay or shale for-mations and silt aquitards [2][3][4], the electrochemical surface effect caused by the water-clay interaction is considered to be the main cause [5,6]. The postlinear flow occurs in the coarse-grained and fractured medium at high velocity; the onset of turbulent flow leads to the inertial force over a viscous force which makes Darcy's law invalid [7,8]. There are two types of functions to explain this nonlinear relationship. One is the Forchheimer equation in which the hydraulic gradient is a second-order polynomial function of the velocity; the first-order term represents the effect of viscous force and the second-order term represents the effect of inertial force. The advantage of this second-order polynomial function is that it has a definite physical meaning and has been validated theoretically [9][10][11]. The other is the Izbash equation in which the hydraulic gradient is a power function of the velocity; this equation is a fully empirical equation based on a large amount of experimental data [12][13][14][15]. Although it lacks a physical meaning, it has been widely used for reasons of mathematical convenience. Both equations can describe non-Darcian to a flow well, and one of the two functions may be favoured in some special cases.
The interlayer staggered zones are special fractures, which have been caused by the localization of shear within relatively narrow zones in fault rocks [16]. They are mainly filled by low-permeability fillings (mud and debris) and the high-permeability fillings (gravel and fractured surrounding rock), and the internal structure is loose, with fractured rocks and high content clay in it [17]. As a result, the interlayer staggered zone has a great influence on the stability of engineering works, such as dam breaks, water inrush in underground excavations, seepage deformation, and leakage in dam foundations [18]. The complex structure of the interlayer staggered zone makes the non-Darcian flow more likely to occur. Moreover, the permeability characteristics of the interlayer staggered zone in non-Darcian flow are still poorly understood.
In situ tests are important ways to estimate the aquifer parameters. There are several categories of in situ hydraulic tests, including constant rate injection or withdrawal tests, slug tests, constant head tests, and recovery tests after constant injection or withdrawal [19]. Except for recovery tests, the non-Darcian flow has been observed in other kinds of the test from many studies. The non-Darcian flow behaviour in a constant rate withdrawal test has been widely studied; some analytical and numerical solutions have also been derived. Sen [20][21][22] firstly used the Boltzmann transform to investigate the non-Darcian flow in a constant pumping test. Dejam et al. [4] provided an exhaustive study of pre-Darcy flow and solved highly nonlinear diffusivity equations with the aid of a generalized Boltzmann transformation technique. Wen et al. [6,15,23,24] developed a linearization procedure to deal with the nonlinear term in the governing equation for constant rate pumping tests in different conditions. Liu et al. [25,26] proposed two generalized non-Darcian radial flow analytical models based on Forchheimer's law and Izbash's law by considering the dimension of flow geometry. Furthermore, the numerical solutions for the non-Darcian flow can deal with more complex flow situations and overcome the disadvantage of the linearization procedure where the true drawdown at early and moderate times [11,[27][28][29] is underestimated. Wang et al. [30] established a non-Darcian flow transient model of slug test in a leaky confined aquifer to quantify the influence of non-Darcian behaviour. Moreover, the non-Darcian flow had been found during several constant head hydraulic tests due to the high hydraulic gradients and velocity [31][32][33]. Chen et al. [13,34,35] proposed two approximate analytical models for data interpretation and estimating hydraulic parameters by performing high-pressure packer tests.
To our knowledge, the unsteady non-Darcian flow in the constant packer test has not been fully studied. Previous studies did not consider the nonlinear flow generated by the constant head packer test that results will be out of step with the actual situation. In this study, a novel analytic solution by taking advantage of the observed flow rate and injection pressure during the constant head packer test is used to describe the unsteady non-Darcian flow in a saturated aquifer. The power law-based non-Darcian flow and a linearization approximation proposed by Wen et al. [6] with the Laplace transform is applied to solve the constant head condition. The model will be verified with the Darcian condition, and a field application is used to obtain the non-Darcian flow parameter from the interlayer staggered stagger zone confined aquifer in the Baihetan hydropower station, China. It provides a reference for further evaluating the seepage stability of the underground building.
The general sketch of this research is shown in Figure 1. First, the theoretical model and its solutions are presented. Then, their results are discussed by a case study and the corresponding model analysis are conducted as well, followed by a summary and conclusions.

Methods
An inverted cone of depression is formed during the packer test, as is shown in Figure 2; the injection pressure can reach a constant high value due to the two expanded inflatable packers to prevent the leaking. Meanwhile, the constant injection pressure can convert to the constant head boundary condition that the water head in the radius of the water injection well is a constant value; the flow rate can also be obtained by the injection system. Therefore, the following assumptions for the packer test can be given that (1) the two-packer includes the whole test section. (2) The whole system is hydrostatic before the test and flow lines are horizontal at the initial moment. ( Summary and conclusions Figure 1: The general sketch of this research. 2 Geofluids where sðr, tÞ is the hydraulic head at radial distance r at time t, qðr, tÞ is the specific discharge, S s is the specific storage, and m is the thickness of aquifer (In this case, it is the thickness of the interlayer staggered zone). With the above assumptions, the initial and boundary conditions can be written as where s w is the injection water head of the injection well with the radius r w which can convert by the injection pressure P (s w = P/ρg). In this study, the power law is adopted to describe the non-Dracian flow where n is a constant which is between 1 and 2 and k 1 is the non-Darcian hydraulic conductivity. Wen et al. [36] proposed the quasihydraulic conductivity k′ to reflect the discharge capacity of the medium. In this study, the relationship between k 1 and k ′ can be expressed as k 1 = ðk ′ Þ 1/n , and the advantage of the expression form k 1 is that it has the same dimension as the hydraulic conductivity which has a clearer physical meaning. Only when n is equal to 1, k 1 is the Darcian hydraulic conductivity. Substituting Equation ( The non-linear term ð∂sðr, tÞ/∂rÞ ðn−1Þ/n should be replaced to obtain a further solution. There are two possible ways to replace the ∂s/∂r term that one is using the finite difference method to obtain a numerical solution and the other is simpler by introducing an approximation term as follows: Equation (5) is the linearization approximation that Q t is the flow rate at each time t. It means the hydraulic gradient at any radial face can be approximated as the relationship between the amount of water passed through at the time t and the non-Darcian hydraulic conductivity k 1 . Note that Q t is an observed constant value at each time point because only in this way the further Laplace transform can be performed next. Fortunately, the flow rate and the constant head can be obtained at the same time by the flowmeter in the packer test equipment. Then, Equation (4) can be converted to a linear partial differential equation: The Laplace transform is an effective measure for solving Equation (6), and the Stehfest algorithm is widely adopted for reversing Laplace transform. For instance, Fan et al. [37,38] adopted this method to deduce analytical solutions for pumping from a poroelastic, confined aquifer considering the effects of a finite-thickness skin zone and the wellbore storage. By using the Laplace transform for time t, the solution of the hydraulic head in the Laplace transform is (see the appendix for derivation) Equation (7) cannot be inverted analytically; one possible way is using the Stehfest method [14,39]. For the most accurate and stable solution, the required number of terms N in the Stehfest method is 18. In this study, a MATLAB program is used to simulate the non-Darcian flow by the above method during the packer test.

Field Application
There are two different ways to analyze the non-Darcian flow model, including the dimensionless form and dimensional form. The dimensionless form is widely used in the study of Darcian flow which is a so-called type curve. The advantage of this method is that it can make full use of observational data to avoid measurement errors in data.

Geofluids
However, in the case of non-Darcian flow, the benefits of the type curves are not obvious due to the nonlinear relationship of the gradient and the discharge [36]. The dimensional form can directly show the sensitivity of the non-Darcian flow parameters (n and k). To analyze the non-Darcian flow behaviours under a high water injection pressure, a set of packer tests were conducted in the interlayer staggered zone aquifer located in the Baihetan hydropower station, southwestern China.
3.1. Site Description. The Baihetan hydropower station is located on the downstream reach of the Jinsha River, China (Figure 3(a)). It is the world's second largest hydropower project which includes a 289 m high dam, a navigation facility, and two underground powerhouses. Figure 3(b) shows the distribution of the underground powerhouse on both sides of the bank which includes (I) one main house, (II) one transformer chamber, (III) four tailrace surge chambers, (IV) eight diversion pipes, and (V) four tailrace tunnels in each bank.
Two underground powerhouses are located in the Basalt rock mass, and the Tuff interlayers are developed between each Basalt layer. The interlayer staggered zones (e.g., C 2 , C 3 , and C 4 ) are formed under structural action in these interlayers (Figure 3(c)); they have higher permeability than the surrounding rocks which represent a potential threat to the stability of the underground powerhouse in future, espe-cially the zone C 2 which runs through the whole underground powerhouse in the left bank.
We conducted packer tests to understand the seepage properties for these interlayer staggered zones. One test site is in the target zone C 2 of the left bank, numbered as CZK57. Zone C 2 is developed in the Tuff layer P 2 β 2 4 which is mainly composed of gravel, rock fragments, and a small amount of medium-dense brown-red mud; it will become soft and plastic when exposed to water. In summary, the permeability of the mud filling is weak, while the gravel and rock fillings are strong [17].

Field
Test. The test site CZK57 is in the tailgate ventilation and safety tunnel of the left bank; Figure 4(a) shows the site conditions in which three test boreholes were set up by using water drilling technology with the elevation of the orifice is 636 m. The distribution of these two boreholes is shown in Figure 4(b), the radius of each borehole is r w = 0:0375 m, and the average thickness of the test section is 0.95 m. As is shown in Figure 4(c), two hydraulic inflatable packers are installed in each borehole to isolate the test section to reach a constant head in the test.
One test was conducted in borehole CZK57 with the injection borehole is CZK57-0, and the observation borehole is CZK57-1; the design injection pressure is 0.3 MPa. In this study, the flow rate (Q) can be obtained by the flowmeter,  5 Geofluids and the water injection pressure (P) can be measured by the pressure gage from the grouting data acquisition and processing system (TIANJIN SAIZHI, FEC-GJ3000). Meanwhile, the TD-Diver groundwater datalogger (DI1271) can record the change of water level in the observation boreholes. These observed data are used in the inversion analysis to get the non-Darcian flow characteristics in the interlayer staggered zone.
The observation data is shown in Figure 5

Parameter Determination.
Combine the aforementioned analytical solution and the test data, the genetic algorithm (GA) is applied to estimate three non-Darcian flow hydraulic parameters: S s , k 1 , and n. The GA is a stochastic search algorithm that provides an efficient search technique for highly nonlinear problems [40]. It has been widely applied   6 Geofluids in groundwater problems [41][42][43][44]. In this study, the decision variables were the three parameters mentioned above, and optimal decision variables which minimize the minimizing the root mean square errors (RMSEs) were found by the genetic algorithm. The RMSE can be written as where m is the number of observation points, h i cal is the hydraulic heads calculated by the model at the ith time, and h i obs is the hydraulic heads measured from the observation borehole at the same ith time.
In this study, the size of the initial population is set to 5000, the function tolerance is 1 × 10 −6 , and the probabilities for crossover and mutation operations are specified to 0.8 and 0.2, respectively. Based on fully investigating the packer test and understanding the hydrogeological conditions of the site, the range of three unknown parameters should be estimated firstly. The composition of the test section is gravel, rock fragments, and mud so that the value of S s should be between 1:0 × 10 −6 m -1 and 1:0 × 10 −4 m -1 by referring to summarized results from Kuang et al. [45]. The non-Darcian index n is between 1 and 2 due to the postlinear flow occurs at high velocity, and the non-Darcian flow permeability k 1 is an unknown parameter that gives a wide range between 1:0 × 10 −7 m·s -1 and 1:0 × 10 −4 m·s -1 .
The genetic algorithm results are depicted in Figure 6; the calculated points in each time can fit the correspondent observation points well, in which the minimum value of RMSE is 0.073 m. The calculated curve is smooth, which represents the observed Q has a slight impact on the change of hydraulic head in the observation borehole. The low permeability of the test section leads to a low flow rate, and the turbulent flow under the non-Darcian flow makes the influence of Q even less obvious. In summary, three non-Darcian flow parameters are listed in Table 1 and the the optimal values are n = 1:278, k 1 = 1:613 × 10 −5 m·s -1 , and S s = 9:757 × 10 −5 m -1 , respectively.

Model Validation.
Another test was used to verify the obtained non-Darcian flow parameters in which the injection borehole is CZK57-1, and the observation borehole is CZK57-2, the distance between these two boreholes is 6 m, and the injection pressure is still 0.3 MPa. Figure 7 shows the recorded flow rate Q and water level h of the borehole CZK57-2; the test lasted 88 minutes. It is similar to the previous test that the flow rate Q increased rapidly at the beginning of the test and then progressively decreased to a stable state. The observed water level of CZK57-2 gradually increased and stabilized with the same growth trend as the previous test.
The inversion results of the previous section are selected to calculate the change of hydraulic head in the observation borehole (r = 6 m). Figure 8 shows the calculated value by our model and the observed value of borehole CZK57-2. The calculated value curve is rough due to the fluctuation of the flow rate. Meanwhile, the calculated value has a good fit with the observed value at the beginning and end of the test, and the observed value is larger during the middle stage of the test. Since there are few observations, we selected the calculated value corresponding to the observation point to calculate the RMSE is 0.091 m. It means the obtained parameters are reliable to describe the non-Darcian flow in the interlayer staggered zone C 2 .

The Distribution of Hydraulic Head in the Non-Darcian
Flow Model. According to the results from the genetic algorithm, the distribution of hydraulic head in the section during the test is shown in Figure 9; six specific time points were selected for display purposes. The hydraulic head is large for small r (less than 0.5 m) due to the constant injection water, and it makes a larger hydraulic gradient around the injection borehole. Other areas have lower hydraulic head and increasing slowly; the influence of non-Darcian flow makes the test section have poor permeability. Moreover, in the first five minutes, the hydraulic head changed more obviously; then, it increases slowly and trends to maintain a stable state after 30 minutes. The emergence of non-Darcian flow will cause a high water head around the injection borehole.  is commonly used to predict the temporal or spatial hydraulic head or drawdown distribution or to determine aquifer parameters for Darcian flow, which can be written as [46,47].

Results and Discussion
where J 0 and Y 0 are Bessel functions of first and second kinds of order 0, respectively, x is a dummy variable, and S is the storativity which S = mS s and T is the trans-missivity of the aquifer which T = mk 1 . The flow rate of the test well can also be obtained from Equation (9) by Darcy's law, where In the comparison of the above model, the value of n is set to 1 in our model which represents the Darcian flow that occurs in this condition, and other parameters are set as r = 3 m, m = 0:95 m, s w = 30 m, S s = 1 × 10 −4 m -1 , and k 1 = 2 × 10 −5 m/s. Meanwhile, the flow rate of the injection well is the calculation results adopted by Equation (10). The comparison of these two models is shown in Figure 10; our model completely agreed with the linear flow model which has only slight differences between them. It means our linearization approximation non-Darcian flow model causes negligible errors and confirms our model is applicable. Figure 11 depicts the sensitivity of the hydraulic head versus n with n = 1, 1.1, 1.2, 1.3, and 1.4, respectively. The hydraulic head in the radius of 3 m decreases when n increases. At the early time of the test, the hydraulic head has a larger increment with the lower n; then, the change of s is slowing down as the test continues. Moreover, the flow approached a steady state more quickly when n gets a larger value, and the hydraulic head becomes lower as well. A greater n value means a greater deviation from the Darcian flow with a greater degree of turbulence flow. Turbulent flow is the most efficient way of dissipating hydraulic energy comparing with laminar flow; therefore, one will observe the slower growth of the hydraulic head for the turbulent flow case. Moreover, it is more difficult for the flow to reach the steady state under higher intensity turbulence conditions.   Figure 12 with S s = 1 × 10 −5 , 5 × 10 −5 , 1 × 10 −4 , 5 × 10 −4 , and 1 × 10 −3 m -1 , respectively. Hydraulic head slow growth for larger specific storage which is similar to that of the power index is demonstrated. When other conditions remain the same, larger specific storage represents that the aquifer has a stronger capacity to take into storage because of head rise from a physical point of constant head packer test. At early times, a larger increment of the hydraulic head is observed with smaller specific storage because of the poor water storage capacity of the aquifer, and it approaches a steady state faster as well. Note that when S s is larger enough, the hydraulic head is not increased rapidly anymore early, and it continually increased in late time.

Effect of the Non-Darcian Hydraulic Conductivity k 1 .
As is shown in Figure 13, the sensitivity of the hydraulic head versus non-Darcian hydraulic conductivity k 1 = 2 × 10 −7 , 2 × 10 −6 , 2 × 10 −5 , 2 × 10 −4 , and 2 × 10 −3 m/s, respectively. Note that the larger k 1 leads to a larger hydraulic head which can represent the permeability of the test section as well, and the trend of hydraulic head increasing is similar to the above sensitivity analysis. This finding can be understood physically. A larger k 1 value indicates a larger "hydraulic conductivity"; it is easier for the medium to transmit water, and the hydraulic head growth is faster early with a larger value of k 1 . As the water injection reaches stability, the flow will reach a quasisteady state faster with better permeability. Additionally, the same situation occurs when k 1 reaches a lower value that the hydraulic head increases slowly in the early time; it shows that the groundwater flow is harder to spread under turbulent flow and weak permeability.

Global Sensitivity Analysis.
To analyze the sensitivity of the model output to the variation of each input parameter, the Sobol method is adopted for global sensitivity analysis. This method has been widely employed in many hydrogeological studies [48][49][50]. The advantage of this method is that no assumption of linearity or monotonicity is required in the adopted interpretative model [51,52]. The Sobol indices are the ratio of the partial model variance to the total model variance which the first-order indices can be written as

Geofluids
And the second-order indices can be written as where V i and V ij represent the partial variance of the model associated with the ith parameter and the interaction of the i th and jth parameters, respectively. V represents the total variance. S i is the first-order sensitivity effect indices that describe the contribution of one parameter to the total variance of the model output. And S ij represents the influence of the interaction of two input parameters on the model output variance. Therefore, the total sensitivity indices S Ti can be expressed as [50] The range of evaluated parameters n, S s , and k 1 is shown in Table 2, and other parameters in the model are the same as the abovementioned value. The whole test process is used to analyze the change of total sensitivity indices. Figure 14 shows the total sensitivity indices versus time during the whole test process. It can be seen that the power index n shows a big influence on the model output variance. The total sensitivity index of n has a larger value at the beginning of the test; then, it increases slightly and tends to remain constant over time. On the other hand, the other two parameters show lower values. The total sensitivity index of k 1 is initially small, then decreases to a lower value with time gradually and trends to remain constant later. The specific storage S s shows a similar trend as the k 1 ; however, it has a higher value k 1 than at the initial time, then sharply decreases within about 1 min and gradually drops to a smaller value than k 1 . It indicates the specific storage S s has little impact on the model output. This is because, at early times, a sudden constant head applied near the   10 Geofluids borehole; the injection water gradually flows into the aquifer. Meanwhile, the specific storage coefficient has an impact on the aquifer. When the injection time is larger enough, the storage is completed and the flow tends to stabilize; turbulence flow becomes the main influence factor at a late time.

Summary and Conclusions
This non-Darcian flow method is simple and feasible; it can be widely used in field tests to descript the distribution of hydraulic head or obtain the corresponding paraments quickly. Meanwhile, it makes full use of observation data in which the variation of flow rate during the test and the injection pressure as the constant boundary. However, the constant head and precise measurements of the test process gives a high requirement for test equipment. One limitation of this method is that only the constant head value of the first stage data can be used in the step-head packer test. In this study, Izbash's law-based linearization procedure presented by Wen et al. [6] was used to derive a new analytic solution for describing the unsteady non-Darcian flow that occurs for the packer test. This method is used to compute the hydraulic head distribution for the packer test conducted in an interlayer staggered zone aquifer system in the Baihetan powerhouse project, China. Then, a set of test data was analyzed, and a genetic algorithm was applied to estimate non-Darcian flow parameters. Finally, the sensitivity of the hydraulic head to three non-Darcian flow parameters was analyzed such as the power index n, the specific storage S s , and the non-Darcian hydraulic conductivity k 1 , and the global sensitivity is also performed. The main findings can be drawn as follows: (1) The flow rate in the injection borehole and the change of water level in the observation borehole in one injection pressure stage of the packer test were used to estimate the non-Darcian flow parameters of the interlayer staggered zone. The estimated results obtained by genetic algorithm are n = 1:278, k 1 = 1:613 × 10 −5 m·s -1 , and S s = 9:757 × 10 −5 m -1 , respectively. It can be used for further evaluating the seepage stability of the underground building (2) The results of this new method for the special Darcian flow case (n = 1) has a perfect fit with the constant head test model solution derived from Darcy's law. Larger values of the power index n and the specific storage S s result in the hydraulic head increases slower during the whole test time. Meanwhile, increasing hydraulic head approaches a steady state slower by taking a larger value. The trend of non-Darcian hydraulic conductivity k 1 is contrary to the previous analysis (3) Each parameter has its influence on the hydraulic head. The head is the most sensitive to the power index n and is insensitive to the specific storage S s and the non-Darcian hydraulic conductivity k 1