Theoretical Model and Experimental Research on Determining Aquifer Permeability Coefficients by Slug Test under the Influence of Positive Well-Skin Effect

Abstract

In the process of slug test in the field, mud and debris may infiltrate into the aquifer during the formation of the test well, and a poorly permeable annulus area appears around the test well wall. This is called the positive well-skin effect, which can have a significant impact on the results of slug tests due to the characteristics of a slug test. Therefore, to carry out the research on the determination of aquifer permeability coefficients by slug test under the influence of positive well-skin effect, a mathematical model of a slug test considering the positive well-skin effect was constructed. The Laplace transform method and the AWG algorithm were used to solve for a series of standard curves, and specific steps were given for calculating aquifer permeability coefficients by slug test under the influence of positive well-skin effect. Meanwhile, an indoor test model with positive well-skin effect was constructed, then several sets of slug tests and pumping tests with different excitation strengths were carried out. Based on this study, the presence of positive well-skin layer had a great impact on the aquifer permeability coefficient determined from both the pumping test and Kipp model of slug test, resulting in smaller results. The HWS model overcomes the influence of positive well-skin effect. This will further improve the theory of slug test, increase the accuracy of slug test in determining the permeability coefficient, and promote the application scope of slug test. It is therefore possible to discern the well-skin factor and determine the aquifer permeability coefficient by fitting the slug test data to a standard curve.

1. Introduction

As an efficient technique for testing the permeability coefficients of aquifers, the slug test is characterized by short test cycle, simple equipment, and high accuracy [1]. Hvorslev [2], Cooper [3], Bouwer and Rice [4], Kipp [5], McElwee [6], Huang [7], Zhou [8], Zhao [9,10,11], and Liang [12] all proposed theoretical models and engineering applications for the slug test in determining hydrogeological parameters of aquifers, none of which have considered the well-skin effect.

However, in the actual application of the slug test, it is likely that a positive well-skin effect appears due to the inevitable influence of factors such as mud wall protection during the drilling. Based on the mechanism characteristics of slug test, this positive well-skin effect greatly affects the accuracy of the slug test and thus has become a technical bottleneck limiting the application of slug test to a large extent. Meanwhile, current methods to reduce the influence of well-skin effect on permeability still focus on how to eliminate the presence of well-skin effect, and while the test results are slightly more accurate, they are costly and inefficient. As a result, various research on the influences of well-skin effect on the slug test have been carried out at home and abroad. Back in 1972, Ramey et al. [13] were the first to propose a theoretical model for slug test that took into account the well-skin effect. They assumed that the well-skin layer was extremely thin, with boundaries extending to infinity, and that the water storage effect of the well was taken into consideration. A series of standard curves were obtained to determine the permeability coefficient of the aquifer and well-skin effect based on the relationship between the dimensionless water storage coefficient and factors related to the well-skin effect. The model considers a more comprehensive situation, but the similarity of the standard curve shapes leads to a high degree of uncertainty in determining the relevant parameters. As a result, this fitting method may lead to inaccurate estimates of aquifer permeability coefficients. Faust [14] investigated the effect of well-skin of finite thickness on the results of the slug test using a numerical model and a simple analytical solution. The results show that the data obtained from the slug test are more representative of the permeability coefficient of the well-skin layer rather than that of the aquifer. The reason is that when well-skin of finite thickness is existent, the curve of the slug test results is shifted along the horizontal axis, making the estimation of the permeability coefficient unreliable. Moench [15] summarized the existing theoretical models of the slug test and pointed out the relationship between the two assumptions of infinitely thin and finitely thick well-skins; the solution for a finite-thickness well-skin is equal to the solution for an infinitely thin well-skin with the storage effect of the well-skin neglected. Sageev [16] assumed well-skin layer existed with extremely thin thickness. Considering the storage effect of the well bore, it was found that when solving for aquifer permeability coefficients with the matching line method, the measured curves should be analysed in three periods: early, middle, and late. The early period mainly reflects the permeability coefficient of well-skin layer, the middle period curve is suitable to the case where there is no well-skin, and the late period curve can be used when well-skins are positive or when the well-skin effect is not considered. Song [17] and Zeng [18] analysed the slug test results under the influence of well-skin effect, divided the slug test process into initial and intermediate stages, and developed a two-stage analysis method to estimate the well-skin factor $S_w$ and aquifer permeability coefficient $K_a$. In a slug test conducted by Rovey II [19] and others on 31 test wells of varying development strength in porous alluvium, a link was found between the well-skin effect and the development strength of the test wells, with the test results being small due to the presence of a positive well-skin effect rather than wellbore damage. The results were also analysed with the means of conventional Bouwer–Rice method and the deconvolution method, which quantifies and removes some degree of well-skin effect. Barrash et al. [20] founded the cause of positive well-skin effect and the evidence of head loss through observation and analysis results of a large number of tests in the BHRS research base: compaction, residue invasion, and blockage of screen. Yeh et al. [21] proposed a simulation method to estimate the thickness of well-skin layer and other hydraulic parameters at the same time and used the annealing (SA) algorithm to estimate five parameters: hydraulic conductivity, specific storage and thickness of the well-skin layer, hydraulic conductivity, and specific storage of the aquifer. Chen [22] first compared the slug test with the conventional pumping test and found that the slug test has very obvious advantages and that it causes less groundwater disturbance, but the radius of influence of slug test is basically within 10 m and only represents a small area of stratigraphic properties near the test well, and the slug test has high requirements on the quality of test well structure and is susceptible to the well-skin effect. To address this drawback, he concluded that the test well structure needs to be pre-designed and adequately washed before the test. Ju [23] carried out pumping test and slug test in the field, and the test results showed that the pumping test involved a large area, while the slug test involved a small area, and the $K$ value of the pumping test was 1.5 times higher than that of the slug test. The well-skin effect was the main reason that leads to the difference between the two test results. Cardiff et al. [24] considered that due to slug tests’ high sensitivity to low-conductivity skins, it was possible that aquifer $K$ estimates obtained may simply represent depth-variability in a low-conductivity skin. By comparing $K$ profiles obtained at closely-spaced wells, the volume being interrogated by the slug test response extended beyond the wellbore skin and represented true aquifer $K$ variability. Sahin [25] proposed that the influence of the well-skin effect must be considered when calculating aquifer permeability parameters in slug tests and that slug test results that took into account the influence of the well-skin effect were more reliable and accurate. Weeks [26] proposed that the well-skin effect had a significant impact on the test results of the slug test, affecting parameters related to the permeability parameters and the water storage coefficient. Liu [27] used numerical approaches to propose a new three-dimensional slug test model for inhomogeneous aquifers in which inertia effects associated with linear friction and well-skin effect are considered. Morozov [28] evaluated aquifer anisotropic permeability and well-skin effect based on slug test data in incomplete wells.

In previous studies, most of the studies on slug test under positive well-skin effect used numerical simulation or parameter inversion of a large number of test data to estimate the aquifer permeability coefficient. There was no good method to accurately calculate the permeability coefficient of well-skin layer, nor was there a special large-scale physical test model to carry out the research of slug test under the positive well-skin effect to determine the aquifer permeability coefficient, so the mechanism research of the well-skin effect on the slug test is lacking. Therefore, it is necessary to study the hydraulic conductivity between the aquifer and the well-skin from the perspective of the mechanism of the well-skin effect. In this study, the mathematical model of the slug test under the positive well-skin effect is firstly proposed; the standard curve of the model was obtained by solving with Laplace transform and AWG (Azari–Wooden–Gaver) algorithm, and the corresponding specific calculation steps of the model based on the solution process are given. Meanwhile, the physical test model was constructed indoors, based on data from multiple sets of slug tests and pumping tests with different excitation strengths, the reliability and accuracy of the proposed theory and method for determining aquifer permeability coefficients by slug test under positive well-skin effects (HHU-Well-Skin model, whose abbreviation is HWS model) was verified by comparing the calculation results with those of the pumping tests and Kipp model method. This will further improve the theory of the slug test, increase the accuracy of the slug test in determining the permeability coefficient, and promote the scope of application of the slug test.

2. Theoretical Model Construction for Slug Test under Positive Well-Skin Effect (HSW Model)

2.1. Basic Condition Assumptions

(1) The aquifer is a confined aquifer of constant thickness and horizontal occurrence, with a homogeneous and isotropic medium; (2) groundwater flowing in the aquifer is radial; (3) there is a well-skin layer of finite thickness, and the distance from the centre of the well to the outer boundary of the well-skin layer is $r_s$; (4) head friction losses are negligible; and (5) the test well is a complete well.

2.2. Establishment and Solution of Model

The continuity equation was established in the confined aquifer and finite thickness well-skin layer; the water balance equation was established in the contact surface between the confined aquifer and the finite thickness well-skin layer and between the finite thickness well-skin layer and the well. The theoretical model of the slug test under the positive well-skin effect was thus constructed, and the aquifer structural model of the slug test under the positive well-skin effect is shown in Figure 1.

(1) Establishment of the theoretical model

In the confined complete well system, the governing equation of groundwater radial flow in a well-skin layer and a confined aquifer is:

$$\frac{{\partial }^2{h}_1}{\partial r^2}+\frac{1}{r}\frac{\partial h_1}{\partial r}=\frac{S_1}{K_1}\frac{\partial h_1}{\partial t}, r_w\le r\le r_s$$

(1)

$$\frac{{\partial }^2{h}_2}{\partial r^2}+\frac{1}{r}\frac{\partial h_2}{\partial r}=\frac{S_2}{K_2}\frac{\partial h_2}{\partial t},r_s\le r\le \infty$$

(2)

A water balance equation is established in the contact surface between the confined aquifer and the well-skin layer of finite thickness:

$$K_1\frac{\partial h_1\left(r_s,t\right)}{\partial r}=K_2\frac{\partial h_2\left(r_s,t\right)}{\partial r},t>0$$

(3)

A water balance equation is established in the contact surface between the well-skin layer of finite thickness and the well:

$$\pi r_c^2\frac{dH\left(t\right)}{dt}=2\pi r_w{K}_{1}B\frac{\partial h_1\left(r_w,t\right)}{\partial r},t>0$$

(4)

Initial conditions:

$$h_1\left(r,0\right)=h_2\left(r,0\right)=0,r>r_w$$

(5)

$$H\left(0\right)=H_0$$

(6)

Boundary conditions:

$$h_2\left(\infty ,t\right)=0$$

(7)

$$h_1\left(r_s,t\right)=h_2\left(r_s,t\right),t>0$$

(8)

$$H\left(t\right)=h_1\left(r_w,t\right)$$

(9)

Among them, subscripts 1 and 2 represent well-skin layer and confined aquifer parameters, respectively; $h$ is the value of head change $\left(m\right); r$ is the distance from a point within the well-skin layer or aquifer to the centre of the well $\left(m\right)$; $t$ is time $\left(s\right)$; $K$ is the permeability coefficient $\left(m·d^{-1}\right)$; $S$ is the water storage coefficient; $r_w$ is the radius of the screen pipe $\left(m\right)$; $r_c$ is the radius of the well casing ($m$); $r_s$ is the distance from the centre of the well to the outer boundary of the well-skin layer $\left(m\right)$; $B$ is the thickness of the confined aquifer $\left(m\right);$ and $H\left(t\right)$ is the change of well water level at time $t$ relative to the initial water level $\left(m\right)$.

(2) Solution of the theoretical model

First, the equations and the definite solution conditions were transformed into dimensionless models by dimensionless factors and parametric variables, and then, in turn, the Laplace transform was performed. The dimensionless transformation equation is:

$$h_{1D}=\frac{h_1}{H_0},h_{2D}=\frac{h_2}{H_0},\rho =\frac{r}{r_w},{\rho }_s=\frac{r_s}{r_w},{\rho }_w=1,\tau =\frac{K_{2}t}{S_2{r}_w^2},\alpha =\frac{K_2}{K_1},\beta =\sqrt{\frac{S_1{K}_2}{S_2{K}_1}},\gamma =\frac{r_c^2}{2S_{2}B{r}_w^2}$$

(10)

The results are as follows after Laplace transformations of time in Equations (1) and (2):

$$\frac{{\partial }^2\overline{h_{1D}}}{\partial {\rho }^2}+\frac{1}{\rho }\frac{\partial \overline{h_{1D}}}{\partial \rho }-{\beta }^{2}p\overline{h_{1D}}=0, 1\le \rho \le {\rho }_s$$

(11)

$$\frac{{\partial }^2\overline{h_{2D}}}{\partial {\rho }^2}+\frac{1}{\rho }\frac{\partial \overline{h_{2D}}}{\partial \rho }-p\overline{h_{2D}}=0, {\rho }_s\le \rho \le \infty$$

(12)

Then, Equations (3) and (4) are, respectively, converted to:

$$\frac{\partial \overline{h_{1D}}}{\partial \rho }\left({\rho }_s,p\right)=\mathsf{\alpha }\frac{\partial \overline{h_{2D}}}{\partial \rho }\left({\rho }_s,p\right), p>0$$

(13)

$$\alpha \gamma \left(p\overline{H_D}-1\right)=\frac{\partial \overline{h_{1D}}\left(1,p\right)}{\partial \rho }$$

(14)

The conditional Equations of definite solution (7)–(9) are, respectively, converted to:

$$\overline{h_{2D}}\left(\infty ,p\right)=0,p>0$$

(15)

$$\overline{h_{1D}}\left({\rho }_s,p\right)=\overline{h_{2D}}\left({\rho }_s,p\right), p>0$$

(16)

According to the modified Bessel equation of order 0,

$$\frac{d^{2}y}{d{x}^2}+\frac{1}{x}\frac{dy}{dx}-{\beta }^{2}y=0$$

(17)

The general solution of Equation (17) is:

$$y=C_1{K}_0\left(\beta x\right)+C_2{I}_0\left(\beta x\right)$$

(18)

Among them, $I_0\left(x\right)$ and $K_0\left(x\right)$are the first and second type of modified Basel functions, respectively; $C_1$and$C_{2 }$are two constants determined by the boundary conditions. Therefore, according to the general solution of Equation (17), the solutions of Equations (11) and (12) can be obtained as:

$$\overline{h_{1D}}\left(\rho ,p\right)=C_1{K}_0\left(\beta \rho \sqrt{p}\right)+C_2{I}_0\left(\beta \rho \sqrt{p}\right)$$

(19)

$$\overline{h_{2D}}\left(\rho ,p\right)=C_3{K}_0\left(\rho \sqrt{p}\right)+C_4{I}_0\left(\rho \sqrt{p}\right)$$

(20)

Among them, $C_1$and $C_2$as well as $C_3$and $C_4$are constants determined by boundary conditions in well-skin layers and aquifer, respectively.

According to Equation (15),

$$\underset{\rho \to \infty }{\lim }\overline{h_{2D}}\left(\rho ,p\right)=0$$

(21)

The first type of modified Bessel function $I_0\left(x\right)$ is an exponentially growing function with the property that when $x\to \infty$, $I_0\left(x\right)\to \infty$, so $C_4=0$ in Equation (21).

$$\overline{h_{2D}}\left(\rho ,p\right)=C_3{K}_0\left(\rho \sqrt{p}\right)$$

(22)

According to Equation (13),

$$C_1{K}_0\left(\beta {\rho }_s\sqrt{p}\right)+C_2{I}_0 \left(\beta {\rho }_s\sqrt{p}\right)=C_3{K}_0\left({\rho }_s\sqrt{p}\right)$$

(23)

For

$$\frac{\partial \overline{h_{1D}}}{\partial \rho }\left({\rho }_s,p\right)=\frac{\partial \overline{h_{1D}}}{\partial \rho }\left(\rho ,p\right){|}_{\rho ={\rho }_s}=C_1\beta \sqrt{p}\frac{\partial K_0\left(\beta {\rho }_s\sqrt{p}\right)}{\partial \left(\beta {\rho }_s\sqrt{p}\right)}+C_2\beta \sqrt{p}\frac{\partial I_0\left(\beta {\rho }_s\sqrt{p}\right)}{\partial \left(\beta {\rho }_s\sqrt{p}\right)}$$

(24)

$$\alpha \frac{\partial \overline{h_{2D}}}{\partial \rho }\left({\rho }_s,p\right)=\alpha \frac{\partial \overline{h_{2D}}}{\partial \rho }\left(\rho ,p\right){|}_{\rho ={\rho }_s}=\alpha C_3\sqrt{p}\frac{\partial K_0\left({\rho }_s\sqrt{p}\right)}{\partial \left({\rho }_s\sqrt{p}\right)}$$

(25)

By the properties of the modified Bessel equation of order 0,

$$\frac{d}{dx}{K}_0\left(x\right)=-K_1\left(x\right)$$

(26)

$$\frac{d}{dx}{I}_0\left(x\right)=I_1\left(x\right)$$

(27)

Therefore, Equations (25) and (26) are:

$$C_1\beta \sqrt{p}\frac{\partial K_0\left(\beta {\rho }_s\sqrt{p}\right)}{\partial \left(\beta {\rho }_s\sqrt{p}\right)}+C_2\beta \sqrt{p}\frac{\partial I_0\left(\beta {\rho }_s\sqrt{p}\right)}{\partial \left(\beta {\rho }_s\sqrt{p}\right) }=-C_1\mathsf{\beta }\sqrt{p}{K}_1\left(\beta {\rho }_s\sqrt{p}\right)+C_2\mathsf{\beta }\sqrt{p} I_1\left(\beta {\rho }_s\sqrt{p}\right)$$

(28)

$$\alpha C_3\sqrt{p}\frac{\partial K_0\left({\rho }_s\sqrt{p}\right)}{\partial \left({\rho }_s\sqrt{p}\right)}=-\alpha C_3\sqrt{p}{K}_1\left({\rho }_s\sqrt{p}\right)$$

(29)

Therefore, according to the Equations (16), (25), (26), (28), and (29), the following conclusions can be obtained:

$$-C_1\beta \sqrt{p}{K}_1\left(\beta {\rho }_s\sqrt{p}\right)+C_2\beta \sqrt{p}{I}_1\left(\beta {\rho }_s\sqrt{p}\right)=-\alpha \sqrt{p}{C}_3{K}_1\left({\rho }_s\sqrt{p}\right)$$

(30)

Simplification of Equation (30):

$$C_1\beta K_1\left(\beta {\rho }_s\sqrt{p}\right)-C_2\beta I_1\left(\beta {\rho }_s\sqrt{p}\right)=\alpha C_3{K}_1\left({\rho }_s\sqrt{p}\right)$$

(31)

For

$$p\overline{H_D}\left(\rho \right)-1=p\overline{h_{1D}}\left(1,p\right)-1=p\left[C_1{K}_0\left(\beta \sqrt{p}\right)+C_2{I}_0\left(\beta \sqrt{p}\right)\right]-1$$

(32)

According to Equations (25) and (29),

$$\frac{\partial \overline{h_{1D}}}{\partial \rho }\left(1,p\right)=\frac{\partial \overline{h_{1D}}}{\partial \rho }\left(\rho ,p\right){|}_{\rho =1}=-C_1\beta \sqrt{p}{K}_1\left(\beta \sqrt{p}\right)+C_2\beta \sqrt{p}{I}_1\left(\beta \sqrt{p}\right)$$

(33)

Associating Equations (32) and (33) according to Equation (14),

$$\alpha \gamma \left[p\left(C_1{K}_0\left(\beta \sqrt{p}\right)+C_2{I}_0\left(\beta \sqrt{p}\right)\right)-1\right]=-C_1\beta \sqrt{p}{K}_1\left(\beta \sqrt{p}\right)+C_2\beta \sqrt{p}{I}_1\left(\beta \sqrt{p}\right)$$

(34)

According to Equation (34),

$$C_1=\frac{\alpha \gamma +C_2\left[\beta I_1\sqrt{p}\left(\beta \sqrt{p}\right)-\alpha \gamma p{I}_0\left(\beta \sqrt{p}\right)\right]}{\alpha \gamma p{K}_0\left(\beta \sqrt{p}\right)+\beta \sqrt{p}{K}_1\left(\beta \sqrt{p}\right)}$$

(35)

Equation (24) $\times \mathsf{\alpha }{K}_1\left({\rho }_s\sqrt{p}\right)$and Equation (31) $\times K_0\left({\rho }_s\sqrt{p}\right),$ respectively:

$$C_1\alpha K_0\left(\beta {\rho }_s\sqrt{p}\right)K_1\left({\rho }_s\sqrt{p}\right)+C_2\alpha I_0\left(\beta {\rho }_s\sqrt{p}\right)K_1\left({\rho }_s\sqrt{p}\right)=C_3\alpha K_0\left({\rho }_s\sqrt{p}\right)K_1\left({\rho }_s\sqrt{p}\right)$$

(36)

$$C_1\beta K_1\left(\beta {\rho }_s\sqrt{p}\right)K_0\left({\rho }_s\sqrt{p}\right)-C_2\beta I_1\left(\beta {\rho }_s\sqrt{p}\right)K_0\left({\rho }_s\sqrt{p}\right)=C_3\alpha K_1\left({\rho }_s\sqrt{p}\right)K_0\left({\rho }_s\sqrt{p}\right)$$

(37)

The solution of combining Equations (36) and (37) is:

$$C_1=\frac{-C_2\left[\alpha I_0\left(\beta {\rho }_s\sqrt{p}\right)K_1\left({\rho }_s\sqrt{p}\right)-\beta I_1\left(\beta {\rho }_s\sqrt{p}\right)K_0\left({\rho }_s\sqrt{p}\right)\right]}{\alpha K_0\left(\beta {\rho }_s\sqrt{p}\right)K_1\left({\rho }_s\sqrt{p}\right)-\beta K_1\left(\beta {\rho }_s\sqrt{p}\right)K_0\left({\rho }_s\sqrt{p}\right)}$$

(38)

From Equations (35) and (38),

$$\begin{array}{c}\left[\alpha I_0\left(\beta {\rho }_s\sqrt{p}\right)K_1\left({\rho }_s\sqrt{p}\right)+\beta I_1\left(\beta {\rho }_s\sqrt{p}\right)K_0\left({\rho }_s\sqrt{p}\right)\right]\times \left[\alpha \gamma p{K}_0\left(\beta \sqrt{p}\right)+\beta \sqrt{p}{K}_1\left(\beta \sqrt{p}\right)\right]C_2\\ -\left[\beta \sqrt{p}{I}_1\left(\beta \sqrt{p}\right)-\alpha \gamma p{I}_0\left(\beta \sqrt{p}\right)\right]\times \left[\beta K_1\left(\beta {\rho }_s\sqrt{p}\right)K_0\left({\rho }_s\sqrt{p}\right)-\alpha K_0\left(\beta {\rho }_s\sqrt{p}\right)K_1\left({\rho }_s\sqrt{p}\right)\right]C_2\\ =\left[\beta K_1\left(\beta {\rho }_s\sqrt{p}\right)K_0\left({\rho }_s\sqrt{p}\right)-\alpha K_0\left(\beta {\rho }_s\sqrt{p}\right)K_1\left({\rho }_s\sqrt{p}\right)\right]\times \alpha \gamma \end{array}$$

(39)

Set:

$$a_1=\alpha \gamma p{I}_0\left(\beta \sqrt{p}\right)-\beta \sqrt{p}{I}_1\left(\beta \sqrt{p}\right)$$

(40)

$$a_2=\alpha \gamma p{K}_0\left(\beta \sqrt{p}\right)+\beta \sqrt{p}{K}_1\left(\beta \sqrt{p}\right)$$

(41)

$$c_1=\alpha I_0\left(\beta {\rho }_s\sqrt{p}\right)K_1\left({\rho }_s\sqrt{p}\right)+\beta I_1\left(\beta {\rho }_s\sqrt{p}\right)K_0\left({\rho }_s\sqrt{p}\right)$$

(42)

$$c_2=\beta K_1\left(\beta {\rho }_s\sqrt{p}\right)K_0\left({\rho }_s\sqrt{p}\right)-\alpha K_0\left(\beta {\rho }_s\sqrt{p}\right)K_1\left({\rho }_s\sqrt{p}\right)$$

(43)

Therefore, according to Equation (39),

$$C_1=\frac{c_1\alpha \gamma }{a_1{c}_2+a_2{c}_1}$$

(44)

According to Equation (38),

$$C_2=\frac{c_2\alpha \gamma }{a_1{c}_2+a_2{c}_1}$$

(45)

According to Equation (32),

$$C_3=\frac{c_1{K}_0\left(\beta {\rho }_s\sqrt{p}\right)+c_2{I}_0\left(\beta {\rho }_s\sqrt{p}\right)}{K_0\left({\rho }_s\sqrt{p}\right)}$$

(46)

Substituting $C_1$, $C_2,$ and $C_3$ into Equations (20) and (21), respectively, the solution of Laplace’s equation is obtained as:

$$\overline{h_{1D}}\left(\rho ,p\right)=\frac{\alpha \gamma \left[c_1{K}_0\left(\beta \rho \sqrt{p}\right)+c_2{I}_0\left(\beta \rho \sqrt{p}\right)\right]}{a_1{c}_2+a_2{c}_1}$$

(47)

$$\overline{h_{2D}}\left(\rho ,p\right)=\frac{\alpha \gamma \left[c_1{K}_0\left(\beta {\rho }_s\sqrt{p}\right)+c_2{I}_0\left(\beta {\rho }_s\sqrt{p}\right)\right]K_0\left(\rho \sqrt{p}\right)}{\left(a_1{c}_2+a_2{c}_1\right)K_0\left({\rho }_s\sqrt{p}\right)}$$

(48)

The final solution of the Laplace transform domain is obtained according to Equation (17) as:

$$\overline{H_D}\left(p\right)=\frac{\alpha \gamma \left[c_1{K}_0\left(\beta \sqrt{p}\right)+c_2{I}_0\left(\beta \sqrt{p}\right)\right]}{a_1{c}_2+a_2{c}_1}$$

(49)

Among them, $\overline{H_D}\left(p\right)$ represents the well water level change in dimensionless form. The rest of the symbols remain the same meaning as above.

Equation (49), which is the final solution derived in the Laplace transform domain, therefore also requires an inverse Laplace transform of Equation (49). However, when using Stehfest’s algorithm [29] for the inversion of the Laplacian space image function, although the analytical solution of the model can be obtained relatively quickly, the number of summation terms $N$ ($N$ is an even number and $N\le 16$) in the algorithm is affected by the Laplacian space solution and the time $t$ and other parameters; especially when the curve of the solution changes steeply, numerical dispersion and oscillation will occur [30]. In recent years, more accurate algorithms such as Crump and AWG have emerged in the study of Laplace numerical inversion methods, with the AWG algorithm being an improvement on the Stehfest algorithm, which maintains the advantages of the original algorithm in terms of simplicity and ease of implementation while ensuring accuracy in solving for steep changes in the curve of the real space solution.

$$g\left(t\right)=\frac{\ln 2}{\tau }{\displaystyle \sum }_{i=1}^N{V}_{i}G\left(\frac{\ln 2}{\tau }i\right)$$

(50)

Among them, $g\left(t\right)$ is the image primitive function in real space, $F\left(s\right)$ is the image function in Laplace space; $N$ is the number of summation terms, $V_i={\left(-1\right)}^{\frac{N}{2}+i}\sum _{k=\left[\frac{i+1}{2}\right]}^{Min\left(i,\frac{N}{2}\right)}\frac{k^{\frac{N}{2}}\left(2k+1\right)!}{\left(\frac{N}{2}-k+1\right)!k!\left(k+1\right)!\left(i-k+1\right)!\left(2k-i+1\right)!}$.

Substitute Equation (49) into Equation (50) for inverse Laplace transform to obtain

$$H_D\left(\tau \right)=\frac{\ln 2}{\tau }{\displaystyle \sum }_{i=1}^N{V}_i\overline{H_D}\left(\frac{\ln 2}{\tau }i\right)$$

(51)

Equation (51) is the analytic solution of the model, where the operation of the parameters $a_1$, $a_2$, $c_1$, and $c_2$ involves the first and second type of modified Bessel functions $I_0\left(x\right)$, $I_1\left(x\right)$ and $K_0\left(x\right)$, $K_1\left(x\right)$, all four of which are in the form of series and are usually simplified according to their properties for the convenience of the operation. Considering that the omission of the remainder term will have an impact on the model calculation results, the Bessel functions besseli(0,x), besseli(1,x) and besselk(0,x), besselk(1,x) are called directly when MATLAB programming is performed, without omitting the remainder term. In addition, to ensure that the solution curve is smooth and realistic, the number of summation terms $N$ in the solution process is taken as 2.

The standard curve data corresponding to different well-skin factors α were obtained by using MATLAB programming to run through Equation (51), and then, the $H_D-\tau$ standard curve under the influence of positive well-skin effect can be obtained by plotting the standard curve data as shown in Figure 2. Under the slug test conditions, the measured curve of ${H^{\prime }}_D-t$ can be plotted according to the slug test data so that the measured curve is aligned with the vertical axis of the standard curve at the same horizontal height position. Then, the measured curve is shifted along the horizontal axis to make the measured curve fit the standard curve, and the corresponding $\alpha$ and $\gamma$ values of the fitted standard curve are recorded so that the α value of the positive well-skin factor of the well can be identified according to the fitted standard curve. At the same time, any matching point is selected, and the corresponding coordinate values $\left[\tau \right]$ and $\left[t\right]$ of the standard curve and the measured curve are noted and substituted into Equations (52) and (53) to obtain:

Aquifer permeability coefficient:

$$K_2=\frac{\left[\tau \right]}{\left[t\right]}\frac{r_c^2}{2\gamma B}$$

(52)

Positive well-skin layer permeability coefficient:

$$K_1=\frac{K_2}{\alpha }$$

(53)

3. Indoor Model Tests

3.1. Test Platform Construction

(1) Aquifer System

Based on the assumptions of the radial flow and the principle of axisymmetry in the groundwater well flow problem, a semi-circular model in the axisymmetric model was utilized for the test platform to save test space. The overall design of the model is a 2.2 m high, 2 m radius semi-cylindrical sand tank with a test well in the centre of the semi-cylindrical section. The test well has an internal diameter of 0.08 m and an external diameter of 0.10 m. Overflow outlets are uniformly arranged at the upper 1.5 m of the model bottom plate on the side of the semi-cylindrical arc to ensure that the boundary is a constant water head boundary during the test. The aquifer system simulates a confined aquifer with a thickness of 0.8 m. The aquifer medium is medium and fine sand with a particle size of 0.25 mm–0.50 mm. Both the top and bottom plates are provided with waterproof layers. A layer of sandbags is laid on the top of the top plate to apply pressure to simulate a confined aquifer. In addition, a gap is left at the circular boundary as a constant water head recharge boundary. As the research focuses on the influence of well-skin on this slug test, the setting of the well-skin layer is very important and two different thicknesses of well-skin layers were used to investigate the influence of different thicknesses of well-skin layer on the results of the slug test. The well-skin layer is located within 20 cm around the test well and is separated from the aquifer using a semi-cylindrical frame, which is divided into two semi-cylindrical zones of 0–10 cm and 10–20 cm. The positive well-skin is filled with a 1:1 mixture of silty clay and silty fine sand. When the well-skin is not set, the well-skin layer is filled with the same particle size of 0.25 mm–0.50 mm sand as the aquifer. The planar graph of the test platform model is shown in Figure 3, where Figure 3a shows the overall planar graph, and Figure 3b shows the local well-skin layer planar graph. The east–west section of the indoor test platform is shown in Figure 4a, and the north–south section is shown in Figure 4b. A physical picture of the indoor test platform is shown in Figure 5.

(2) Observation System

In order to ensure the reliability of the observation data, the observation system uses both piezometric pipes and automatic water pressure acquisition sensor. The main function of this system is to observe the change of water level data in real time during the test. The test well S0 is set in the centre of the semi-cylindrical section as part of both the piezometric pipe observation device and the automatic water pressure acquisition sensor device, and the piezometric pipes and automatic water pressure acquisition sensors are distributed in two radial directions of the semi-circular base plate.

The piezometric pipe observation device consists of the test well S0 and piezometric pipes C1#–C9#, which are, respectively, 0.05 m, 0.10 m, 0.15 m, 0.20 m, 0.25 m, 0.50 m, 1.00 m, 1.50 m, and 2.00 m away from the test well. The automatic water pressure acquisition sensor device consists of the test well S0 and observation holes S1#–S8#, which are, respectively, 0.05 m, 0.10 m, 0.15 m, 0.25 m, 0.50 m, 1.00 m, 1.50 m, and 2.00 m away from the test well. An automatic water pressure acquisition sensor was put into the test well S0 from the wellhead to record the water level change of the test well. The overall setup of the observation system is shown in Figure 3.

3.2. Test Programme

(1) Test Programme Design

The slug test is a method to determine the hydrogeological parameters of an aquifer by means of the change law of well water level with time caused by the change of small amount of water in the test well. For this indoor slug test, three different volumes of water lifting columns were used as excitation devices to meet the requirements of different excitation strengths: small (500 cm³), medium (1000 cm³), and large (1500 cm³) to meet the instantaneous characteristic of the water volume changes in the test well. In addition, to eliminate the possible influence of factors such as positive well-skin layer thickness and excitation strength, the indoor test programme was divided into three different well-skin conditions: no well-skin, 10 cm thick positive well-skin layer, and 20 cm thick positive well-skin layer. Three different excitation strengths for the water lifting slug test and three different pumping flow rates for large, medium, and small flow rates pumping test were conducted under different well-skin conditions. To avoid accidental errors, three parallel sets of tests were conducted simultaneously for each test. The pumping test and slug test numbers for each group are shown in

Table 1. Table of test numbers.

Pumping Flow Rate (Q/m3·d−1)	Pumping Test Method
	No Well-Skin	10 cm Thick Positive Well-Skin Layer	20 cm Thick Positive Well-Skin Layer
		${\mathit{\rho }}_{\mathit{s}}=3.75$	${\mathit{\rho }}_{\mathit{s}}=6.25$
Small flow rate (12.4416 m3·d−1)	N-Ps1	P10-Ps1	P20-Ps1
	N-Ps2	P10-Ps2	P20-Ps2
	N-Ps3	P10-Ps3	P20-Ps3
Medium flow rate (17.4528 m3·d−1)	N-Pm1	P10-Pm1	P20-Pm1
	N-Pm2	P10-Pm2	P20-Pm2
	N-Pm3	P10-Pm3	P20-Pm3
Large flow rate (23.3280 m3·d−1)	N-Pl1	P10-Pl1	P20-Pl1
	N-Pl2	P10-Pl2	P20-Pl2
	N-Pl3	P10-Pl3	P20-Pl3
Excitation strength		Slug test method
	No well-skin	10 cm thick positive well-skin layer	20 cm thick positive well-skin layer
		${\rho }_s$= 3.75	${\rho }_s$= 6.25
Small column (500 cm3)	N-Sls1	P10-Sls1	P20-Sls1
	N-Sls2	P10-Sls2	P20-Sls2
	N-Sls3	P10-Sls3	P20-Sls3
Medium column (1000 cm3)	N-Slm1	P10-Slm1	P20-Slm1
	N-Slm2	P10-Slm2	P20-Slm2
	N-Slm3	P10-Slm3	P20-Slm3
Large column (1500 cm3)	N-Sll1	P10-Sll1	P20-Sll1
	N-Sll2	P10-Sll2	P20-Sll2
	N-Sll3	P10-Sll3	P20-Sll3

.

(2) Test Steps and Methods

Before the start of the test, water was injected into the aquifer from the side of the model arc to saturate it layer by layer, with the air in the sand body drained outward. Water injection was stopped when the water level in the model completely immersed the aquifer sand body. The sand body was considered reaching saturation when the readings of each piezometric pipe and automatic water pressure acquisition sensor were basically the same and remained unchanged for a long time. Then, water was injected continuously until the water level reached the overflow outlets of the model, and there was water overflow. The pumping test and slug test were carried out after the water level in the model was always at the overflow outlets position, and the readings of each piezometric pipe and automatic water pressure acquisition sensor were stable and remained unchanged for a long time.

4. Test Results and Analysis

4.1. Pumping Test Method

The results of the pumping tests were processed and calculated using the Dupuit equation method [31] and the Thiem equation method [31]. The average value of the results of the two calculation methods for each group of pumping tests without well-skin model conditions was used as the reference value to verify the reliability and accuracy of the proposed theory and method (HWS model) for determining the aquifer permeability coefficient by slug tests under the positive well-skin effect. The results of the permeability coefficients calculated by the pumping test method without well-skin model conditions are shown in

Table 2. Calculated results of pumping tests to determine permeability coefficients under the no well-skin model conditions.

Test Group	Dupuit Equation K Avg./cm·s−1	Theim Equation K Avg./cm·s−1	K Avg./cm·s−1
N-Ps1/2/3	0.1239	0.1197	0.1218
N-Pm1/2/3	0.1144	0.1095	0.1120
N-Pl1/2/3	0.1108	0.1074	0.1091
K Avg./cm·s−1	0.1164	0.1121	0.1143

.

By comparing the mean values of the results calculated by the two methods in Table 2, they are basically consistent, indicating a reasonable test design. Therefore, the average value of the aquifer permeability coefficient calculated under both methods was used as the reference value of the aquifer permeability coefficient determined by the pumping tests without well-skin conditions.

After setting up the positive well-skin layer, the pumping test observation data from the positive well-skin layer and the aquifer were calculated using the Theim equation method to obtain the reference value of the permeability coefficient of the positive well-skin layer. The calculation results obtained are shown in

Table 3. Calculation of results permeability coefficients determined by pumping tests under positive well-skin model conditions.

Test Group	Dupuit Equation K Avg./cm·s−1	Theim Equation K Avg./cm·s−1
	Based on Test Well Test Data	Based on Aquifer Observation Data	Based on Positive Well-Skin Layer Observation Data
P10-Ps1/2/3	0.0274	0.0960	0.0250
P10-Pm1/2/3	0.0222	0.0775	0.0280
P10-Pl1/2/3	0.0188	0.0880	0.0298
K Avg./cm·s−1	0.0228	0.0872	0.0276
P20-Ps1/2/3	0.0182	0.0549	0.0166
P20-Pm1/2/3	0.0170	0.0627	0.0176
P20-Pl1/2/3	0.0130	0.0664	0.0174
K Avg./cm·s−1	0.0161	0.0613	0.0172

.

As shown in Table 3, the aquifer permeability coefficient determined from pumping tests became significantly smaller after the installation of the positive well-skin layer. In the case of the 20 cm thick well-skin, the permeability coefficient of the aquifer determined from pumping tests was even smaller due to the smaller permeability coefficient of the well-skin layer, indicating that the presence of the positive well-skin layer had a great impact on the aquifer permeability coefficient determined from pumping tests.

However, the permeability coefficient of the positive well-skin layer can be accurately obtained by using the observation data in the positive well-skin layer. Therefore, the calculation results of the observation data in the positive well-skin layer using Theim’s equation method show that the reference values of the permeability coefficient in the positive well-skin layer as ${\rho }_s$ = 3.75 and ${\rho }_s$ = 6.25 are 0.0276 cm·s⁻¹ and 0.0172 cm·s⁻¹, respectively.

4.2. Slug Test Method

The slug test data were analysed and calculated with the Kipp model and the proposed theoretical model for slug test under positive well-skin effect (HWS model), respectively.

(1) Kipp model analysis method [5] (without considering the well-skin effect model)

The Kipp model uses dimensionless factors and parameter variables to convert the differential equation of water flow oscillations within a single well into dimensionless form and obtain a set of standard curves. The test data were processed, and the permeability coefficient values were solved using the matching method, and the partial Kipp model matching fitting diagram without well-skin model conditions are shown in Figure 6, and the final calculation results are shown in

Table 4. Calculation results of permeability coefficients of the Kipp model in the slug test.

None Well-Skin		${\mathit{\rho }}_{\mathit{s}}=3.75$		${\mathit{\rho }}_{\mathit{s}}=6.25$
Test Group	K Avg./cm·s−1	Test Group	K Avg./cm·s−1	Test Group	K Avg./cm·s−1
N-Sls1/2/3	0.1583	P10-Sls1/2/3	0.0704	P20-Sls1/2/3	0.0452
N-Slm1/2/3	0.1456	P10-Slm1/2/3	0.0594	P20-Slm1/2/3	0.0554
N-Sll1/2/3	0.1267	P10-Sll1/2/3	0.0528	P20-Sll1/2/3	0.0528
K Avg./cm·s−1	0.1435		0.0609		0.0511

.

From Table 4, the average values of the permeability coefficients calculated from the slug tests with different excitation strengths under the no well-skin model conditions are in general agreement with the reference values of the permeability coefficients calculated from the conventional pumping tests. This indicates that accurate results can be obtained from the slug test data using the Kipp model under the no well-skin model conditions. However, the calculated values based on the Kipp model under positive well-skin model conditions are between the aquifer permeability coefficient and the positive well-skin layer permeability coefficient, which represents the average value of the positive well-skin layer and aquifer permeability coefficient under the influence of positive well-skin. This implies that the slug test is influenced by positive well-skin, and the Kipp model is no longer applicable to the slug test under positive well-skin model conditions.

(2) Analytical approach to the theoretical model of the slug test under positive well-skin effect (HWS model)

According to the well-skin layer model conditions of the indoor test, standard curves were plotted for no well-skin (${\rho }_s$ = 1, $\beta$ = 1, $\gamma$ = 1000), 10 cm thick positive well-skin layer (${\rho }_s$ = 3.75, $\beta$ = 5, $\gamma$ = 1000), and 20 cm thick positive well-skin layer (${\rho }_s$ = 6.25, $\beta$ = 5, $\gamma$ = 1000), respectively. We next took the test data of $H_D$ and time $t$ in the slug test under different well-skin model conditions and drew a single logarithmic relationship curve, fitted it with the standard curve (Figure 7, Figure 8 and Figure 9), chose any matching point, and recorded the coordinate values of $\left[\tau \right]$ and $\left[t\right]$ corresponding to the standard curve and the measured curve, and then the aquifer permeability coefficient $K_2$ could be obtained. Then, based on the value of $\alpha$, the well-skin layer permeability coefficient $K_1$ could be obtained.

The HWS model was used to analyze the test data of each group. The calculation results of the aquifer and well-skin layer permeability coefficients were obtained, respectively, and shown in

Table 5. Calculation results of permeability coefficient determined by slug test based on HWS model under different well-skin effects.

Test Number	$\mathit{\alpha }$	Aquifer Permeability Coefficient K2 Avg./cm·s−1
N-Sls1/2/3	1	0.1330
N-Slm1/2/3	1	0.1110
N-Sll1/2/3	1	0.1430
K Avg./cm·s−1		0.1290
Test number	$\alpha$	Aquifer permeability coefficient K2 Avg./cm·s−1	Positive well-skin layer permeability coefficient K1 Avg./cm·s−1
P10-Sls1/2/3	5	0.1250	0.0250
P10-Slm1/2/3	5	0.1250	0.0250
P10-Sll1/2/3	5	0.1100	0.0220
K Avg./cm·s−1		0.1200	0.0240
P20-Sls1/2/3	8	0.1100	0.0140
P20-Slm1/2/3	9	0.1250	0.0140
P20-Sll1/2/3	8	0.1000	0.0130
K Avg./cm·s−1		0.1120	0.0140

.

4.3. Comparison and Analysis of Calculation Results

Based on the calculation results in Table 2, Table 3, Table 4 and Table 5, the permeability coefficients of aquifer and well-skin layers in different well-skin conditions obtained by different calculation methods were compared and analysed as shown in Figure 10, Figure 11 and Figure 12.

The average value of the Dupuit equation method, Theim equation method, Kipp model, and HWS model was used as the reference value of aquifer permeability coefficient (0.1252 cm·s⁻¹) under the no well-skin model conditions. The calculation results of different test methods show that the calculation results of Kipp model of the slug test are generally larger than the reference value, and the calculation results of HWS model are more similar to the reference value; that is, the HWS model can obtain accurate aquifer permeability coefficient values based on the slug test data under the no well-skin model. Meanwhile, it can be seen from Figure 10 that the calculation results based on the slug test are more stable under the same excitation strength.

Comparing the results of the different analytical methods for ${\rho }_s$ = 3.75 (10 cm thick positive well-skin layer), both the pumping test and the slug test Kipp model analytical methods are affected by the positive well-skin effect in the calculation of the aquifer permeability coefficient, resulting in small results. The HWS model overcomes the influence of the positive well-skin effect. The average value (0.1200 cm·s⁻¹) of the HWS model is in good agreement with the reference value (0.1252 cm·s⁻¹) of permeability coefficient. Meanwhile, the average value (0.0240 cm·s⁻¹) of the permeability coefficient of the positive well-skin layer calculated by the HWS model is consistent with the average value (0.0276 cm·s⁻¹) of the permeability coefficient of the positive well-skin layer calculated using the Theim equation method based on the observation test data within the positive well-skin layer. This shows that the permeability coefficient obtained from the HWS model for ${\rho }_s$ = 3.75 is accurate under positive well-skin effect.

Comparing the results obtained by different analytical method for ${\rho }_s$ = 6.25 (20 cm thick positive well-skin layer), it can be seen that as the thickness of the well-skin layer becomes larger, the positive well-skin effect has a greater impact on the calculation results of the pumping test. Compared to the pumping test and the Kipp model of the slug test, the average value (0.1120 cm·s⁻¹) of the aquifer permeability coefficient calculated by the HWS model for ${\rho }_s$ = 6.25 is essentially the same as the reference value of aquifer permeability coefficient (0.1252 cm·s⁻¹) under the no well-skin model conditions. The calculation results of the permeability coefficient of the positive well-skin layer for ${\rho }_s$ = 6.25 of both the pumping test Theim equation method and the HWS model are smaller than those for ${\rho }_s$ = 3.75, probably because the thickening of the positive well-skin layer makes the filling medium more compact and leads to a smaller permeability coefficient value. The average value (0.0140 cm·s⁻¹) calculated by the HWS model for the positive well-skin layer of ${\rho }_s$ = 6.25 is similar to the average value (0.0172 cm·s⁻¹) of the permeability coefficient of the positive well-skin layer calculated using the Theim equation method based on the observation test data within the positive well-skin layer, indicating that the HWS model still maintains a high accuracy in calculating the permeability coefficient values for the aquifer and the positive well-skin layer at ${\rho }_s$ = 6.25.

Comparing the calculations results of aquifer permeability coefficient under different well-skin conditions, the Dupuit equation method of pumping test is more influenced by the positive well-skin effect, as the results obtained when using the test well for calculations represent the average permeability coefficient of the whole aquifer (including positive well-skin layer). Under the positive well-skin effect, the Theim equation method can obtain accurate values of the permeability coefficient of the positive well-skin layer based on the observation test data of the pumping test within the positive well-skin layer, but the calculated permeability value of the aquifer is greatly affected by the positive well-skin effect, and the calculation result is smaller. Compared to the Kipp model, which is only applicable to the no well-skin model conditions, the HWS model is applicable to various positive well-skin models and can accurately obtain the permeability coefficients of aquifer under various positive well-skin model conditions. The permeability coefficient of the positive well-skin layer calculated by the HWS model is similar to the value calculated by the Theim equation based on the observations test data of the pumping test within the positive well-skin layer. This implies that HWS model is accurate in the application of determining the permeability coefficient of aquifer using slug test data, and the HWS model can calculate the permeability coefficient value of the positive well-skin layer when the positive well-skin effect exists. Compared with the traditional pumping test method and Kipp model method of slug test, HWS model is more accurate in determining the aquifer permeability coefficient.

The proposed slug test theoretical model under the influence of positive well-skin effect was verified and compared in indoor finite scale model. Therefore, the application of the proposed theoretical model and method in the field test will be the focus of future research work. In the field application, it is necessary to consider the influence of scale effect, different excitation strength, and different excitation modes (water injection mode, water lift mode, and air pressure mode, etc.) on the test results.

5. Conclusions

• The theoretical model of the slug test considering the positive well-skin effect was solved by using the Laplace transform method and the AWG algorithm. Furthermore, multiple sets of standard curves under different well-skin conditions were plotted, and the specific parameter calculation methods and steps of the HWS model were proposed.
• The presence of the positive well-skin layer had a great impact on the aquifer permeability coefficient determined from both pumping test and slug test, resulting in smaller results. The Kipp model is no longer applicable to the slug test under the influence of positive well-skin effect. The HWS model can overcome the influence of the positive well-skin effect, and the HWS model is applicable to the various positive well-skin models and the no well-skin model. When the positive well-skin effect exists, the HWS model can calculate not only the permeability coefficient of aquifers but also the permeability coefficient of the positive well-skin layer.
• Multiple groups of standard curves under different well-skin conditions have obvious curve characteristics, and different standard curves have a high degree of discrimination. Therefore, by analysing the curve shape and characteristics of the HWS model, according to the different degree of the positive well-skin effect, it can be judged in the field test whether the positive well-skin effect exists based on the preliminary understanding of formation lithology.