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

Authors: Qi Shen; Zhifang Zhou; Meng Chen; Sijia Li; Zhe Wang; Yabing Li
Journal: Geofluids
 (2021)

---

## Abstract

Groundwater flow in an aquifer has frequently been found to be non-Darcian by performing in situ tests. A novel analytic model is proposed in this study for describing the unsteady non-Darcian flow in a confined aquifer by taking advantage of the observed flow rate and injection pressure during the constant head packer test. A linearization approximation of the Izbash equation is used to approximate the nonlinear term in the governing equation. This analytic model is applied to describe the non-Darcian flow in the interlayer staggered zone at the Baihetan hydropower station, China. The test results inversed by the genetic algorithm show that non-Darcian flow happened during the test under the injection pressure 0.3 MPa with the power indexn is 1.278, non-Darcian hydraulic conductivity k1 is 1.613×10−5 cm/s and the specific storage Ss is 9.757×10−5 m-1, respectively. The sensitivity analysis indicated that when the power index n or the specific storage Ss is larger, and the hydraulic head will increase more slowly and needs longer to stabilize, but the non-Darcian hydraulic conductivity k1 shows the opposite trend. Moreover, the hydraulic head is more sensitive to the power index n compared to other parameters at late times. The findings of this study reveal the non-Darcian flow during the constant head packer test and provide a simple and fast way to estimate parameters for more accurate seepage field simulation.


---

## Body

## 1. 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 low-permeability 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 low-permeability 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 formations and silt aquitards [2–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–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–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–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–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–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 Figure1. 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.Figure 1
The general sketch of this research.


## 2. Methods

An inverted cone of depression is formed during the packer test, as is shown in Figure2; 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. (3) The test aquifer is horizontal, homogeneous, isotropic, uniformly thick, and of laterally infinite extent. Then, the continuity equation can be expressed as [6]
(1)∂qr,t∂r+qr,tr=Ss∂sr,t∂t,where sr,t is the hydraulic head at radial distance r at time t, qr,t is the specific discharge, Ss 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
(2)sr,0=0,s∞,t=0,srw,t=sw,where sw is the injection water head of the injection well with the radius rw which can convert by the injection pressure P (sw=P/ρg). In this study, the power law is adopted to describe the non-Dracian flow
(3)qr,t=k1∂sr,t∂r1/n,where n is a constant which is between 1 and 2 and k1 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 k1 and k′ can be expressed as k1=k′1/n, and the advantage of the expression form k1 is that it has the same dimension as the hydraulic conductivity which has a clearer physical meaning. Only when n is equal to 1, k1 is the Darcian hydraulic conductivity. Substituting Equation (3) to Equation (1) yields
(4)∂2sr,t∂r2+nr∂sr,t∂r=Ssnk1∂sr,t∂rn−1/n∂sr,t∂t.Figure 2
The schematic diagrams of the problem.The non-linear term∂sr,t/∂rn−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:
(5)∂sr,t∂r=qr,tk1n≈Qt2πrmk1n.Equation (5) is the linearization approximation that Qt 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 k1. Note that Qt 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:
(6)∂2sr,t∂r2+nr∂sr,t∂r=Ssnk1nQt2πmn−1r1−n∂sr,t∂t.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)
(7)s¯r,p=swrwn−1/2r−n−1/2pK1−n/3−n2Ap2‐n/3−nrw3−n/2K1−n3−n2Ap2‐n3−nr3−n/2.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.


## 3. 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. 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 (Figure3(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.Figure 3
The location of the Baihetan hydropower station and the test sites. (a) Geographical location map of Baihetan hydropower station in China. (b) The layout of hydropower station and the distribution of the test sites. (c) The spatial distribution of main hydrogeological structures (sections A-A’).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.,C2, C3, and C4) 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, especially the zone C2 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 zoneC2 of the left bank, numbered as CZK57. Zone C2 is developed in the Tuff layer P2β24 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].


### 3.2. Field Test

The test site CZK57 is in the tailgate ventilation and safety tunnel of the left bank; Figure4(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 rw=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.Figure 4
The distribution and sectional view of the testing site of borehole CZK57: (a) photo of the test site, (b) the distribution of three boreholes, and (c) the profile of the three boreholes (all dimensions are in units of m).
(a)(b)(c)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, 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 Figure5 under the injection pressure is 0.3 MPa including the change of flow rate and the observation water level. For this test, the distance between the injection borehole CZK57-0 and the observation borehole CZK57-1 is 3 m. The initial water level of the observation borehole is 10.79 m; then, it increased continuously under the constant injection pressure. This injection process lasted 81 minutes, and there is a big fluctuation of flow rate at the beginning of the test, then gradually decreased and stabilized to an average of 52.36 L/min (8.73×10−4 m3/s) 10 minutes after the start of the test. It suddenly dropped at 40 minutes and then stabilized at 50.63 L/min (8.44×10−4 m3/s).Figure 5
Flow rate of injection borehole and water level of observation borehole CZK57-1.


### 3.3. 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:Ss, k1, 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 in groundwater problems [41–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
(8)RMSE=1m∑i=1mhical−hiobs2,where m is the number of observation points, hical is the hydraulic heads calculated by the model at the ith time, and hiobs 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 is1×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 Ss 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 k1 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 Figure6; 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, k1=1.613×10−5 m·s-1, and Ss=9.757×10−5 m-1, respectively.Figure 6
Measured observation borehole hydraulic headh (black points) and calculated hydraulic head by the non-Darcian flow model (red line) of the first test.Table 1
The genetic algorithm inversion of non-Darcian flow parameter results.
MethodNon-Darcian flow parametersRMSE (m)nk1 (m·s-1)Ss (m-1)Genetic algorithm1.2781.613×10−59.757×10−50.073


### 3.4. 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.Figure7 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.Figure 7
Flow rate of injection borehole and water level of observation borehole CZK57-2.The inversion results of the previous section are selected to calculate the change of hydraulic head in the observation borehole (r=6 m).Figure8 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 C2.Figure 8
Measured observation borehole hydraulic headh and calculated hydraulic head by the non-Darcian flow model of the second test for validating obtained non-Darcian flow parameters.


### 3.5. 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 Figure9; 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.Figure 9
The distribution of hydraulic head during the whole packer test in the test section by the non-Darcian flow model.


## 3.1. Site Description

The Baihetan hydropower station is located on the downstream reach of the Jinsha River, China (Figure3(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.Figure 3
The location of the Baihetan hydropower station and the test sites. (a) Geographical location map of Baihetan hydropower station in China. (b) The layout of hydropower station and the distribution of the test sites. (c) The spatial distribution of main hydrogeological structures (sections A-A’).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.,C2, C3, and C4) 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, especially the zone C2 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 zoneC2 of the left bank, numbered as CZK57. Zone C2 is developed in the Tuff layer P2β24 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].


## 3.2. Field Test

The test site CZK57 is in the tailgate ventilation and safety tunnel of the left bank; Figure4(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 rw=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.Figure 4
The distribution and sectional view of the testing site of borehole CZK57: (a) photo of the test site, (b) the distribution of three boreholes, and (c) the profile of the three boreholes (all dimensions are in units of m).
(a)(b)(c)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, 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 Figure5 under the injection pressure is 0.3 MPa including the change of flow rate and the observation water level. For this test, the distance between the injection borehole CZK57-0 and the observation borehole CZK57-1 is 3 m. The initial water level of the observation borehole is 10.79 m; then, it increased continuously under the constant injection pressure. This injection process lasted 81 minutes, and there is a big fluctuation of flow rate at the beginning of the test, then gradually decreased and stabilized to an average of 52.36 L/min (8.73×10−4 m3/s) 10 minutes after the start of the test. It suddenly dropped at 40 minutes and then stabilized at 50.63 L/min (8.44×10−4 m3/s).Figure 5
Flow rate of injection borehole and water level of observation borehole CZK57-1.


## 3.3. 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:Ss, k1, 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 in groundwater problems [41–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
(8)RMSE=1m∑i=1mhical−hiobs2,where m is the number of observation points, hical is the hydraulic heads calculated by the model at the ith time, and hiobs 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 is1×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 Ss 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 k1 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 Figure6; 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, k1=1.613×10−5 m·s-1, and Ss=9.757×10−5 m-1, respectively.Figure 6
Measured observation borehole hydraulic headh (black points) and calculated hydraulic head by the non-Darcian flow model (red line) of the first test.Table 1
The genetic algorithm inversion of non-Darcian flow parameter results.
MethodNon-Darcian flow parametersRMSE (m)nk1 (m·s-1)Ss (m-1)Genetic algorithm1.2781.613×10−59.757×10−50.073


## 3.4. 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.Figure7 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.Figure 7
Flow rate of injection borehole and water level of observation borehole CZK57-2.The inversion results of the previous section are selected to calculate the change of hydraulic head in the observation borehole (r=6 m).Figure8 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 C2.Figure 8
Measured observation borehole hydraulic headh and calculated hydraulic head by the non-Darcian flow model of the second test for validating obtained non-Darcian flow parameters.


## 3.5. 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 Figure9; 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.Figure 9
The distribution of hydraulic head during the whole packer test in the test section by the non-Darcian flow model.


## 4. Results and Discussion

### 4.1. Comparison of the Non-Darcian Flow Model with the Constant Head Test Model

The constant head test model 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].
(9)sr,t=sw1−2π∫0∞exp−TSx2tY0rxJ0rwx−J0rxY0rwxJ02rwx+Y02rwxdxx,where J0 and Y0 are Bessel functions of first and second kinds of order 0, respectively, x is a dummy variable, and S is the storativity which S=mSs and T is the transmissivity of the aquifer which T=mk1. The flow rate of the test well can also be obtained from Equation (9) by Darcy’s law, where
(10)Qrw,t=8Tswπ∫0∞exp−TSx2t1J02rwx+Y02rwxdxx.In the comparison of the above model, the value ofn 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, sw=30 m, Ss=1×10−4 m-1, and k1=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 10
Comparison of the curve for the non-Darcian flow model and constant head test model.


### 4.2. Effect of the Power IndexN

Figure11 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 11
Hydraulic head in the aquifer (r=3 m) with m=0.95 m, sw=30 m, Ss=1×10−4 m-1, and k1=2×10−5 m/s for the values of power index n=1, 1.1, 1.2, 1.3, and 1.4, respectivity.


### 4.3. Effect of the Specific StorageSs

The sensitivity of the hydraulic head versus specific storageSs is shown in Figure 12 with Ss=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 Ss is larger enough, the hydraulic head is not increased rapidly anymore early, and it continually increased in late time.Figure 12
Hydraulic head in the aquifer (r=3 m) with m=0.95 m, sw=30 m, n=1.2, and k1=2×10−5 m/s for the values of power index Ss=1×10−5, 5×10−5, 1×10−4, 5×10−4, and 1×10−3 m-1, respectivity.


### 4.4. Effect of the Non-Darcian Hydraulic Conductivityk1

As is shown in Figure13, the sensitivity of the hydraulic head versus non-Darcian hydraulic conductivity k1=2×10−7, 2×10−6, 2×10−5, 2×10−4, and 2×10−3 m/s, respectively. Note that the larger k1 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 k1 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 k1. As the water injection reaches stability, the flow will reach a quasisteady state faster with better permeability. Additionally, the same situation occurs when k1 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.Figure 13
Hydraulic head in the aquifer (r=3 m) with m=0.95 m, sw=30 m, n=1.2, and Ss=1×10−4 m-1 for the values of power index k1=2×10−7, 2×10−6, 2×10−5, 2×10−4, and 2×10−3 m/s, respectivity.


### 4.5. 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–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
(11)Si=ViV.And the second-order indices can be written as(12)Sij=VijV,where Vi and Vij represent the partial variance of the model associated with the ith parameter and the interaction of the ith and jth parameters, respectively. V represents the total variance. Si is the first-order sensitivity effect indices that describe the contribution of one parameter to the total variance of the model output. And Sij represents the influence of the interaction of two input parameters on the model output variance. Therefore, the total sensitivity indices STi can be expressed as [50]
(13)STi=Si+∑j≠iSij+⋯.The range of evaluated parametersn, Ss, and k1 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.Table 2
Ranges of parameters used in the global sensitivity analysis.
Uncertain parameterRangePower indexn1-2Specific storageSs (m-1)1×10−7-1×10−4Non-Darcian hydraulic conductivityk1 (m/s)1×10−6-1×10−4Figure14 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 k1 is initially small, then decreases to a lower value with time gradually and trends to remain constant later. The specific storage Ss shows a similar trend as the k1; however, it has a higher value k1 than at the initial time, then sharply decreases within about 1 min and gradually drops to a smaller value than k1. It indicates the specific storage Ss has little impact on the model output. This is because, at early times, a sudden constant head applied near the 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.Figure 14
The total sensitivity indices for three uncertain parameters as a function of time.


## 4.1. Comparison of the Non-Darcian Flow Model with the Constant Head Test Model

The constant head test model 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].
(9)sr,t=sw1−2π∫0∞exp−TSx2tY0rxJ0rwx−J0rxY0rwxJ02rwx+Y02rwxdxx,where J0 and Y0 are Bessel functions of first and second kinds of order 0, respectively, x is a dummy variable, and S is the storativity which S=mSs and T is the transmissivity of the aquifer which T=mk1. The flow rate of the test well can also be obtained from Equation (9) by Darcy’s law, where
(10)Qrw,t=8Tswπ∫0∞exp−TSx2t1J02rwx+Y02rwxdxx.In the comparison of the above model, the value ofn 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, sw=30 m, Ss=1×10−4 m-1, and k1=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 10
Comparison of the curve for the non-Darcian flow model and constant head test model.


## 4.2. Effect of the Power IndexN

Figure11 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 11
Hydraulic head in the aquifer (r=3 m) with m=0.95 m, sw=30 m, Ss=1×10−4 m-1, and k1=2×10−5 m/s for the values of power index n=1, 1.1, 1.2, 1.3, and 1.4, respectivity.


## 4.3. Effect of the Specific StorageSs

The sensitivity of the hydraulic head versus specific storageSs is shown in Figure 12 with Ss=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 Ss is larger enough, the hydraulic head is not increased rapidly anymore early, and it continually increased in late time.Figure 12
Hydraulic head in the aquifer (r=3 m) with m=0.95 m, sw=30 m, n=1.2, and k1=2×10−5 m/s for the values of power index Ss=1×10−5, 5×10−5, 1×10−4, 5×10−4, and 1×10−3 m-1, respectivity.


## 4.4. Effect of the Non-Darcian Hydraulic Conductivityk1

As is shown in Figure13, the sensitivity of the hydraulic head versus non-Darcian hydraulic conductivity k1=2×10−7, 2×10−6, 2×10−5, 2×10−4, and 2×10−3 m/s, respectively. Note that the larger k1 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 k1 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 k1. As the water injection reaches stability, the flow will reach a quasisteady state faster with better permeability. Additionally, the same situation occurs when k1 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.Figure 13
Hydraulic head in the aquifer (r=3 m) with m=0.95 m, sw=30 m, n=1.2, and Ss=1×10−4 m-1 for the values of power index k1=2×10−7, 2×10−6, 2×10−5, 2×10−4, and 2×10−3 m/s, respectivity.


## 4.5. 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–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
(11)Si=ViV.And the second-order indices can be written as(12)Sij=VijV,where Vi and Vij represent the partial variance of the model associated with the ith parameter and the interaction of the ith and jth parameters, respectively. V represents the total variance. Si is the first-order sensitivity effect indices that describe the contribution of one parameter to the total variance of the model output. And Sij represents the influence of the interaction of two input parameters on the model output variance. Therefore, the total sensitivity indices STi can be expressed as [50]
(13)STi=Si+∑j≠iSij+⋯.The range of evaluated parametersn, Ss, and k1 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.Table 2
Ranges of parameters used in the global sensitivity analysis.
Uncertain parameterRangePower indexn1-2Specific storageSs (m-1)1×10−7-1×10−4Non-Darcian hydraulic conductivityk1 (m/s)1×10−6-1×10−4Figure14 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 k1 is initially small, then decreases to a lower value with time gradually and trends to remain constant later. The specific storage Ss shows a similar trend as the k1; however, it has a higher value k1 than at the initial time, then sharply decreases within about 1 min and gradually drops to a smaller value than k1. It indicates the specific storage Ss has little impact on the model output. This is because, at early times, a sudden constant head applied near the 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.Figure 14
The total sensitivity indices for three uncertain parameters as a function of time.


## 5. 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 Ss, and the non-Darcian hydraulic conductivity k1, 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 aren=1.278, k1=1.613×10−5 m·s-1, and Ss=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 Ss 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 k1 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 indexn and is insensitive to the specific storage Ss and the non-Darcian hydraulic conductivity k1



---

*Source: 7170769-2021-11-08.xml*