*Article* **Numerical Simulation of Phosphorus Release with Sediment Suspension under Hydrodynamic Condition in Mochou Lake, China**

**Yu Bai 1, Jinhua Gao 2,\* and Tianyi Zhang <sup>3</sup>**


Received: 25 December 2018; Accepted: 18 February 2019; Published: 21 February 2019

**Abstract:** Phosphorus is a major cause of lake eutrophication. Understanding the characteristics regarding the release of phosphorus from sediments under hydrodynamic conditions is critical for the regulation of lake water quality. In this work, the effects of sediment suspension on the release characteristics of phosphorus from sediment were investigated under different hydrodynamic conditions. The experimental results showed that in the experimental process, the phosphorus was at first released quickly into the overlying water but then slowed down. Furthermore, the process of dissolved phosphorus (DP) release under hydrodynamic conditions with and without sediment suspension was simulated using a lattice Boltzmann method. The simulation showed satisfying results.

**Keywords:** dissolved phosphorus; hydrodynamic condition; Lattice Boltzmann method; release characteristics

#### **1. Introduction**

Sewage, industrial discharge, and agricultural runoff contribute to most sources of phosphorus released into shallow lakes, and most of the phosphorus is accumulated in lake sediments, which come back to overlying water under certain conditions [1,2]. Investigations on large and shallow lakes [3,4] indicated that wind waves are critical to sediment resuspension processes. Dynamic behaviors of suspended sediments and wind–wave effects have been studied in past decades and has shown that the suspended sediment has a strong effect on the release of phosphorus [5]. Phosphorus is a major factor that leads to the eutrophication of lakes as it can be absorbed by the vegetation of a lake, and studies have shown that strong wave disturbance can double the phosphorus concentration in the overlying water of lakes [6,7]. Therefore, understanding the release characteristics of the phosphorus from sediments under hydrodynamic conditions is critical for the regulation of lake water quality.

Two methods were proposed to investigate the influence of hydrodynamic forces on the release of phosphorus from sediment. One is to measure the concentration of phosphorus in natural lakes under natural hydrodynamic conditions, and then to analyze the effect of the hydrodynamic force on the phosphorus release from the sediment. The other is to take sediment samples from lakes and then simulate the hydrodynamic conditions as natural conditions using laboratory instruments [8]. In terms of the research on natural lakes, Nilsson et al. measured the phosphorus concentration in the overlying water of sediment under different hydrodynamic conditions in recirculating channel and found that the change of hydrodynamic force could significantly affect the phosphorus concentration in overlying water [9]. Previous research has shown that sediment concentration has a significant

correlation with nutrients [10]. Wang et al. [11] found a positive correlation between total phosphorus (TP) and total suspended sediment (TSS) concentrations. Moreover, numerical methods have been developed to simulate the phosphorus concentrations in lakes and rivers [12]. Some of these methods consider the processes of adsorption–desorption and the release from bed sediment [13].

In recent years, the lattice Boltzmann method (LBM) has been widely applied to the study of diffusion. Wolf-Gladrow [14] gave the basic equations and the parameter selection method of LBM for solving diffusion problems. Jiaung et al. [15] used the LBM to simulate the process of heat conduction, where the thermal diffusivity was controlled by changing the relaxation time. Chen showed that the LBM is also applicable to simulate the diffusion problem in porous media with uniform and heterogeneous porosities, such as lake sediment [16].

In this study, laboratory experiments were conducted to investigate the release characteristics of phosphorus from sediment in shallow lakes. Experiments with and without sediment suspension were designed. The main objectives were as follows: (1) to investigate the release characteristics of the TP and dissolved phosphorus (DP) from the sediment with and without sediment suspension, and (2) to simulate the DP release with and without sediment suspension using LBM.

#### **2. Methodology**

#### *2.1. Laboratory Experiments*

The experiment was divided into two parts: a phosphorus release experiment with sediment suspension and a phosphorus release experiment without sediment suspension. The difference between them was the diameter of sediment, which, depending on its size, would either suspend or not under a hydrodynamic condition.

#### 2.1.1. Experimental Device

The structure of the experimental device is shown in Figure 1a,b. A Plexiglas container with a diameter of 30 cm and a height of 50 cm was mounted with a variable speed motor that ran the propeller to adjust the hydrodynamic conditions in the device.

#### 2.1.2. Experimental Sediment

The sediment was collected from Site 1 and Site 2 (Figure 2), respectively. Due to the previous lake slope protection project, a large amount of large diameter sediment was filled in the shore of Mochou Lake in 2015, the sediment diameter of Site 1 along the shore of Mochou Lake is larger than that of Site 2. The diameter of the sediment at Site 2 (D2) ranged from 0.15 cm to 0.23 cm, with a median diameter of 0.18 cm. The sediment was too coarse to suspend, even at the highest propeller speed (300 rad/min). The sediment diameter (D1) of Site 1 is given in Table 1. The diameter of sediment in Site 1 is small enough to have it suspend in a hydrodynamic condition.



Mochou Lake has a history of 1500 years, with an average water depth of 2.5 m, a lake surface area of 0.3 km2, and a maximum water depth of 4 m. On 21 July 2018, we used a Petersen grab to dig a surface sediment at a depth of 0 cm to 10 cm in Mochou Lake. This lake is frequently influenced by winds, with wind-induced currents transporting dissolved matters (e.g., nutrients). There are two main parts of the sampling device: a sampling grab and a pulling rope. The sampling grab is made of a high-quality alloy material, as shown in Figure 1c. The one-time sampling sediment volume was 1–5 L, and the diameter of the sampling grab was 18 cm.

**Figure 1.** The structure of experimental device and sampling device. (**a**) structure of the experimental device; (**b**) Sketch of the experimental device; (**c**) sampling grab.

**Figure 2.** Sample location and structure of the experimental device.

#### 2.1.3. The Experimental Method

Before the experiment, the sediment samples were fully stirred, and the thickness was controlled at 10 cm. A total of 20 L of deionized water was slowly added to the Plexiglas container, and the container was placed for 24 h prior to the experiment. The speed of the propeller for the various steps was set to 100 rad/min, 200 rad/min, and 300 rad/min, respectively. The experimental parameters are shown in Table 2.

Three 100 g sediment samples were taken from the sediment bed before and after each experiment run. The TP concentration of the sediment was determined by Chinese environmental standard HJ632-2011: Soil Determination of Total Phosphorus by alkali fusion-Mo-Sb anti-spectrophotometry method.

For each experiment run, three 30 mL water samples were extracted from each experimental run for each hour of total experimental time. The DP and TP in the water was determined using the Mo-Sb anti-spectrophotometry method of Yuan [17], the samples for TP determination needed to be digested in advance, and the samples for DP determination needed to be filtered without being digested. In addition, the particulate phosphorus (PP) was estimated as the difference between the TP and the DP (i.e., PP = TP − DP). The TSS concentration was calculated using a drying method. The beaker with 20 mL sampling water was placed in an oven at 115 ◦C till all the water evaporated. The calculated equation was as follows:

$$S = \frac{\mathcal{W}\_1 - \mathcal{W}\_2}{V} \tag{1}$$

where *W*<sup>1</sup> is the weight after drying, *W*<sup>2</sup> is the weight of the beaker, and *V* is the volume of sampling water.


**Table 2.** Experimental parameters.

#### 2.1.4. Selection of Propeller Speeds Based on Field Data

At Site 1, the flow velocity (*uz*), positioned in the water 0.5 m above the bed, was measured using a LS1206B intelligent portable velocimeter which was made by Nanjing Shunlaida Measurement and Control Equipment Co., Ltd. (Nanjing, China). The velocity measurement error was less than 1.5%, and each sampling time was 5 min. The wind speed (*uw*), positioned 0.2 m above the water surface, was measured using an AR866 handheld thermosensitive anemometer made by Suzhou R.B.T Measurement and Control Technology Co., Ltd. (Suzhou, China). The wind speed measurement error was less than 1%, the measurement range was from 0.3 m/s to 30 m/s, and each sampling time was 5 min. The data was measured at 10:00–12:00 am, 2:00–4:00 pm on 22, 24 and 27 July 2018, respectively. Three groups of corresponding velocity and wind speed per hour were measured at Site 1. A statistical analysis of the field data (Figure 3) indicates that the occurrences of the *uz* were less than 0.25 m/s and that the *uw* was less than 8 m/s.

Based on the similarity principle, the rotational speeds are selected so that the Froude number in the laboratory condition (*Frm*) is equal to that in field conditions (*Frp*) [18].

$$Fr\_m = Fr\_p \tag{2}$$

Under laboratory conditions,

$$Fr\_m = \frac{u\_m}{\sqrt{gh\_m}}\tag{3}$$

where *um* is the tangential speed in the laboratory experiment, *hm* is the height above the sediment surface in the laboratory experiment (herein, *hm* = 0.25 m), and *g* is the gravitational acceleration.

Under the field condition

$$Fr\_p = \frac{u\_z}{\sqrt{gh\_z}}\tag{4}$$

where *hz* is the height above the sediment surface in the field (herein, *hz* = 0.5 m) under laboratory conditions.

From the general relation between the tangential and angular speeds, the laboratory experimental rotational speed is

$$
\mu\_{\rm II} = \frac{P\_s}{60} r \tag{5}
$$

where *Ps* is the propeller speed, and *r* is the radius of the stirring rod (herein, *r* = 0.045 m).

This study examined the relationship between *um* and *uw* that corresponds to the dynamic condition caused by the propeller speed to the wind speed. When measuring, the maximum wind speed can be less than 8 m/s, but sometimes the maximum wind speed in Nanjing can be larger than 8 m/s. The field wind speed varied from 4.65 m/s to 12.74 m/s, and *Ps* varied from 0.075 to 0.225 (Table 3). As a result, the propeller speeds were set to 100, 200, and 300 rad/min in the laboratory experiment to simulate the prevailing wind speed of 4.65 m/s to 12.74 m/s in the lake. Huang et al. [18] examined the shear stress caused by lake currents, because they are considered to have similar dynamic effects to the stirring rod on the sediment resuspension, and because the dynamic simulator can produce shear stresses of water currents in the laboratory. The propeller speeds were selected to simulate the prevailing bottom flow velocities of 0 m/s to 0.08 m/s (calculated by Equation (6)), as in Taihu Lake. As a result, the blade stirrer operated at propeller speeds of 0, 100 rad/min, 200 rad/min, 300 rad/min, and 400 rad/min in the laboratory experiment. We used the same method to examine which propeller speeds to select to simulate the prevailing bottom flow velocities of 0 m/s to 0.04 m/s in the lake. As a result, the blade stirrer operated at the propeller speeds of 0, 100 rad/min, 200 rad/min, and 300 rad/min in this laboratory experiment.

The bottom boundary velocity is computed as

$$
\mu\_b = \frac{\kappa u\_z}{\ln\left(\frac{z}{z\_0}\right)}\tag{6}
$$

where *κ* is the von Kármán's constant, and *z*<sup>0</sup> is the bottom boundary roughness, (we assume this based on existing literature [1,8], *z*<sup>0</sup> = 0.02 m).

**Figure 3.** The relation between field flow velocity and wind speed.

**Table 3.** The contrast table of wind velocity and rotational speeds computed.


Note: *um* was calculated by Equation (4); *uw* was calculated by the Equation in Figure 3; *uz* was calculated by Equations (1)–(3).

#### *2.2. Simulation Method*

Generally, the nutrient release in lakes undergoes the following process: the nutrient is released from sediment, the resuspended sediment is desorbed, and then it is diffused in overlying water. For sediment release, it follows the process of diffusion, adsorption, and desorption.

#### 2.2.1. LBM

LBM includes two steps: migration and collision. In the migration step, the particles move in a certain direction and at a certain speed to the adjacent nodes of the grid. The migration process is expressed as:

$$f\_k(\mathbf{x}, y, t + \Delta t) = f\_k(\mathbf{x}, y, t) + f\_k'(\mathbf{x}, y, t) \tag{7}$$

where *fk* is the particle distribution function in terms of a discrete particle along the direction k; *f <sup>k</sup>* is the value of a particle before a migration along the direction k.

Theoretically, the particle collision process is very complex and difficult to solve. Bhatnagar et al. proposed a simple BGK collision operator for a discretized LBM equation, which can be expressed as:

$$f\_k'(\mathbf{x}, y, t) = -\omega f\_k(\mathbf{x}, y, t) + \omega f\_k^{c\eta}(\mathbf{x}, y, t) \tag{8}$$

where *f eq <sup>k</sup>* is the distribution of the equilibrium function along the direction *k*, and *ω* is the relaxation frequency.

Bringing Equation (6) into Equation (5), we have:

$$f\_k(\mathbf{x}, y, t + \Delta t) = f\_k(\mathbf{x}, y, t)[1 - \omega] + \omega f\_k^{eq}(\mathbf{x}, y, t) \tag{9}$$

In LBM, the solution region is divided into many lattices. In this paper, two dimensions with a nine-direction (D2Q4) lattice configuration were used (Figure 4).

$$f\_k^{\text{eq}}(\mathbf{x}, z, t) = w\_k T(\mathbf{x}, z, t) \tag{10}$$

where *fk* is the particle distribution function in terms of a discrete particle along the direction *k*, *f eq k* is the equilibrium distribution function along the direction *k*; w1 = w2 = w3 = w4 = 0.25, and *T* is the nutrients concentration.

**Figure 4.** Nutrient migration process in the experimental device.

#### 2.2.2. Phosphorus Release from Sediment

The phosphorus release in lakes includes the phosphorus release from sediment and the phosphorus diffusion in overlying water. For sediment release, it follows the process of diffusion, adsorption, and desorption. The diffusion rate can be measured by the nutrient concentration gradient

*Water* **2019**, *11*, 370

in the pore water of sediment with Fick's first law, and the adsorption and desorption can be defined by a source item [19].

The phosphorus release from sediment can be expressed as:

$$
\rho \frac{\partial \mathbb{C}}{\partial t} = \wp D\_{zs} \frac{\partial^2 c\_d}{\partial z^2} - \rho\_b \frac{\partial c\_s}{\partial t} \begin{pmatrix} for \ z \ < \ 0 \end{pmatrix} \tag{11}
$$

where *cd* is the DP concentration in water (mg/L), *t* is time (s), *ϕ* is the porosity of sediment, *z* is vertical axis originated (at the sediment-water interface, *z* = 0), *Dzs* is the diffusion coefficient of phosphorus in sediment (m2/s), *ρ<sup>b</sup> ∂cs <sup>∂</sup><sup>t</sup>* is a source term for phosphorus adsorption and desorption by the sediment bed [20], *cs* is the quantity of phosphorus adsorption in the sediment bed (mg/kg) and *ρ<sup>b</sup>* is the density of sediment (kg/m3).

The Lagergren first-order (LFO) equation is commonly used for describing the adsorption and desorption and for their kinetics research, which is expressed as [21]:

$$
\rho\_b \frac{\partial c\_s}{\partial t} = \rho\_b b\_1 (c\_s - c\_{\text{sc}}) \tag{12}
$$

where *b*<sup>1</sup> is the first-order rate constant: absorption efficiency of sediment bed (g/(Ls)), and *cse* is the sediment contamination level (mg/kg).

Yuan et al. [22] assumed that the desorption amount of the sediment bed are equal to the amount added to the solution. Then, they modified the LFO equation as:

$$
\rho\_b \frac{\partial c\_s}{\partial t} = \rho\_b b\_1 (c\_d - c\_{\mathfrak{e}}) \tag{13}
$$

where *ce* is the equilibrium concentration of TP in water (mg/L) [23].

The modified LFO equation only considers the constant hydrodynamic condition where *ce* remains unchanged in a constant condition. Therefore, in an airtight container without phosphorus input, the concentration of phosphorus in overlying water and sediment would be constant values. However, under the action of shear velocity, *ce* will vary with the change of the hydrodynamic condition (or the adsorption rate decreases and *ce* increases with the increase of shear velocity [24]).

Here the second term on the right of Equation (1) is modified as:

$$
\rho\_b \frac{\partial \mathbf{c}\_s}{\partial t} = \rho\_b b\_1 (\mathbf{c}\_d - \mathbf{c}\_c) \tag{14}
$$

*Dzs* can be expressed as:

$$D\_{zs} = \boldsymbol{\varphi}^{m-1} D\_{zm} \tag{15}$$

where *Dzm* is the molecular diffusion coefficient in water (m2/s) and it varies with the targeting solution; *m* = 3 is a constant [25].

The source term can be added in the LBM function [26], and the LBM function can be described as:

$$f\_k(\mathbf{x}, z, t + \Delta t) = f\_k(\mathbf{x}, z, t)(1 - \omega\_p) + \omega\_p f\_k^{eq}(\mathbf{x}, z, t) - \frac{\Delta t w\_k \rho\_b b\_1}{\varrho} (c\_d(\mathbf{x}, z, t) - c\_\varepsilon) \tag{16}$$

where Δ*t* is the time step and *ω<sup>p</sup>* is the relaxation frequency in pore water;

The relationship between the diffusion coefficient and *ω<sup>p</sup>* can be obtained by Chapman-Enskog expansion [26].

$$D\_{zs} = \frac{\Delta x^2}{p \Delta t} \left(\frac{1}{\omega\_p} - 0.5\right) \tag{17}$$

where *p* presents the dimension of the model. (*p* = 1 presents one dimension, *p* = 2 presents two dimensions, and *p* = 3 presents three dimensions).

In this paper, we adopted a two-dimension model, and *ω* can be obtained as follows:

$$
\omega\_p = 1/\left(\frac{2\Delta t D\_{zs}}{\Delta x^2} + 0.5\right) \tag{18}
$$

#### 2.2.3. Nutrient Diffusion in Overlying Water

In overlying water, the formulations can be simply described as a diffusion process, as the phosphorus release from TSS and the biochemical reactions are negligible. The governing equations can be expressed as [27]:

$$\frac{\partial c\_d}{\partial t} = (D\_{zt} + D\_{zm})\frac{\partial^2 c\_d}{\partial z^2} - S\frac{\partial c\_p}{\partial t} \quad (\text{for } z \ge 0) \tag{19}$$

$$\frac{D\_{zt}}{v} = \left(A\frac{zu\_\*}{v}\right)^n\tag{20}$$

where *Dzt* is the turbulent diffusion coefficient (m2/s), *A* is the area of water–sediment interface (m2), *<sup>v</sup>* is the kinematic viscosity of water (m2/s), and *<sup>n</sup>* = 3 is a constant [28]; *<sup>S</sup> <sup>∂</sup>cp <sup>∂</sup><sup>t</sup>* is a source term for phosphorus adsorption and desorption by TSS, *S* is the TSS concentration in the overlying water, and *cp* is the quantity of phosphorus adsorption by TSS in the overlying water.

We assumed that the desorption amount of TSS is equal to the amount added to the solution. Then, *<sup>S</sup> <sup>∂</sup>cp <sup>∂</sup><sup>t</sup>* can be modified as:

$$S\frac{\partial \mathcal{L}\_p}{\partial t} = b\_2 S(\mathbf{c}\_d - \mathbf{c}\_\varepsilon) \tag{21}$$

where *b*<sup>2</sup> is the first-order rate constant, that is, the absorption efficiency of TSS (g/(Ls)).

The LBM function and *ω<sup>w</sup>* can be obtained with the same method as with the sediment:

$$f\_k(\mathbf{x}, z, t + \Delta t) = f\_k(\mathbf{x}, z, t)(1 - \omega\_w) + \omega\_w f\_k^{eq}(\mathbf{x}, z, t) - \Delta t w\_k b\_2 \mathbf{S} \left(\mathbf{c}\_d(\mathbf{x}, z, t) - \mathbf{c}\_t\right) \tag{22}$$

$$
\omega\_w = 1/\left(\frac{2\Delta t(D\_{zt} + D\_{zm})}{\Delta x^2} + 0.5\right) \tag{23}
$$

where *u*<sup>∗</sup> is the shear velocity, and the *u*<sup>∗</sup> generated by the propeller speed is given by Chandler [29]:

$$
\mu\_\* = 8.67 \times 10^{-5} P\_s - 3.27 \times 10^{-3} \tag{24}
$$

#### 2.2.4. Nutrient Desorption from Resuspended Sediment

In this research, the sediment concentration in the overlying water was assumed to be uniform. Here, the dynamic release amount caused by sediment suspension is much larger than the constant release amount caused by diffusion [30], while the vertical velocity component is relatively smaller than the convection velocity component [31]. Therefore, only the effects of sediment re-suspension and sedimentation are considered, and Equation (4) could be simplified into:

$$\frac{dS}{dt} = (S\_r - S\_s)/H \tag{25}$$

$$S\_{\mathcal{T}} = \begin{cases} \begin{array}{c} k(\mathfrak{u}\_{\ast} - \mathfrak{u}\_{\mathcal{E}}) & \mathfrak{u}\_{\ast} > \mathfrak{u}\_{\mathcal{E}} \\ 0 & \mathfrak{u}\_{\ast} < \mathfrak{u}\_{\mathcal{E}} \end{array} \end{cases} \tag{26}$$

$$S\_s = \omega\_d S \tag{27}$$

The settling velocity of sediment particles was obtained by Song et al. [32]:

$$
\omega\_d = \frac{\nu}{D} D\_\*^3 \left( 38.1 + 0.93 D\_\*^{12/7} \right)^{-7/8} \tag{28}
$$

$$D\_\* = \left(\frac{\Delta \mathbf{g}}{\nu^2}\right)^{1/3} D \tag{29}$$

where *<sup>D</sup>* is the particle diameter and *<sup>D</sup>*<sup>∗</sup> is the dimensionless particle diameter. <sup>Δ</sup> <sup>=</sup> *<sup>ρ</sup><sup>s</sup> <sup>ρ</sup>* − 1, where *ρ<sup>s</sup>* and *<sup>ρ</sup>* represent the density of particles and the density of the fluid, respectively. K = 9 × <sup>10</sup>−<sup>3</sup> (kg/m3); *Ss* is the deposition flux; *Sr* is the re-suspension flux; *uc* is the critical velocity of sediment, which can be expressed as:

$$
\mu\_{\varepsilon} = \sqrt{d\_{cr} \Delta \lg D\_{\ast}} \tag{30}
$$

Van Rijn [21] related the dimensionless particle parameter to a critical sediment mobility parameter *(dcr*), which is elaborated in Table 4.


**Table 4.** Variation of *dcr* with *D*∗.

The LBM function for sediment in overlying water can be described as follows,

$$\mathcal{S}\_{(t+\Lambda t)} = \frac{\mathcal{S}\_{r(t)} - \mathcal{S}\_{s(t)}}{H} + \mathcal{S}\_{(t)} \tag{31}$$

#### 2.2.5. Boundary Conditions

By applying LBM, the free surface of overlying water can be defined as the thermal insulating boundary. The nutrient concentration gradient of the free surface is 0, and the container wall is defined as rebound boundary. Please refer to Mohamad [26] for an in-depth explanation of the boundary method.

#### 2.2.6. Model Parameters

The parameters in this model are given in Table 5.


**Table 5.** Parameters in the model.

To conclude the above contents, the sediment and phosphorus transport can be calculated using LBM, as shown in Figure 5.

**Figure 5.** Flow chart of 2-D LBM.

#### **3. Experimental Results**

For Run 4–6, there is very little sediment in the samples, and all of the measured TSS concentration are less than 0.02 g/L, which means that the phosphorus absorbed by TSS in the overlying water can be ignored. Therefore, only the TSS concentration of Run 1–3 are listed in Table 6. The TSS concentration is higher with the higher propeller speed.

The TP concentrations are shown in Table 7; the TP concentrations of Run 1 to Run 6 ranged from 0.059 mg/L to 0.119 mg/L, 0.063 mg/L to 0.191 mg/L, 0.065 mg/L to 0.24 mg/L, 0.031 mg/L to 0.128 mg/L, 0.04 mg/L to 0.143 mg/L, and 0.031 mg/L to 0.151 mg/L, respectively, for the 8 experimental hours. The TP concentration in the water increased quickly in the first four hours, then slowed down in the last four hours. The DP concentrations are shown in Table 8; the DP concentrations of Run 1 to Run 6 ranged from 0.048 mg/L to 0.104 mg/L, 0.034 mg/L to 0.163 mg/L, 0.034 mg/L to 0.207 mg/L, 0.037 mg/L to 0.123 mg/L, 0.039 mg/L to 0.136 mg/L, and 0.03 mg/L to 0.146 mg/L, respectively, for the 8 experimental hours. The DP concentration in the water increased quickly in the first four hours, then slowed down in the last four hours.


**Table 6.** The TSS concentration under different hydrodynamic conditions (g/L).

The variations of TP concentration in the sediment before and after experiment are shown in Figure 6. We can see that the initial TP concentration of Run 1–3 is higher than Run 4–6. After the experiment TP, decreased by 18.5%, 29.7%, 38.8%, 36.7, 46.2%, and 48.9%, from Run 1 to Run 6, respectively.

**Figure 6.** Variations of TP concentration in the sediment before and after the experiment.

#### **4. Simulation Results**

The simulation result of TSS concentration is shown in Figure 7. From Figure 7, we can see that the predicted value fit well with the measured data.

**Figure 7.** Comparison of predicted and measured values of TSS (S is the TSS concentration).

Error analysis was conducted to determine the difference between the predicted and measured data. The average absolute error Δ*<sup>a</sup>* and average relative error Δ*<sup>r</sup>* were calculated by:

$$\overline{\Delta\_d} = \frac{1}{N} \sum\_{1}^{N} abs \left( c\_d \right)\_{computed} - \left( c\_d \right)\_{measured} \right) \tag{32}$$

$$\overline{\Delta\_r} = \left[ \frac{1}{N} \sum\_{1}^{N} abs \left( \frac{(c\_d)\_{\text{computed}} - (c\_d)\_{\text{measured}}}{(c\_d)\_{\text{measured}}} \right) \right] \times 100\% \tag{33}$$

where *N* is the number of measured lines in the cross section for each case, and (*cd*)*computed* and (*cd*)*measured* are the calculated and measured values of the DP concentration.

The DP predicted results are shown in Table 9, along with the measured data. The average relative errors Δ*r* are 10.85%, 8.41%, and 13.01% for simulations 1 to 3, respectively. The average relative errors Δ*r* are 11.67%, 8.84%, and 3.76% for simulations 4 to 6, respectively. Both the suspension model and the non-suspension model show satisfying results. The comparison of results between actual measurement and simulation verified the effectiveness of the LBM.


**Table 9.** Comparison between predicted value and measured value: DP concentration.

#### **5. Conclusions**

The release characteristics of the TP and DP from sediment were experimentally investigated under different hydrodynamic conditions, with and without sediment suspension. By analyzing the measured and predicted data, the main results can be summarized as follows:

In the case of sediment suspension and no-sediment suspension, at the beginning of the experiment, TP and DP were quickly released into the overlying water, and then gradually slowed down. In this experiment, the propeller speeds corresponding to the wind speed in the field condition could improve our result in the field work. LBM was verified to be effective in simulating the process of DP release from sediment under hydrodynamic conditions within experimental conditions, and this technique could provide some theoretical references for the application of the Mochou Lake simulation.

**Author Contributions:** Conceptualization, Y.B. and J.G.; Methodology, Y.B.; Software, Y.B.; Validation, T.Z., Y.B. and J.G.; Formal analysis, T.Z.; Investigation, J.G.; Resources, J.G.; Data curation, J.G.; Writing—original draft preparation, Y.B.; Writing—review and editing, Y.B.

**Funding:** This research received no external funding.

**Acknowledgments:** The authors thank the anonymous reviewers, Associate Editor and Editor for their constructive comments.

**Conflicts of Interest:** The authors declare no conflict of interest.

#### **Notation**


#### **References**


© 2019 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (http://creativecommons.org/licenses/by/4.0/).

## *Article* **The Dynamics of Water Wells Efficiency Reduction and Ageing Process Compensation**

**Krzysztof Polak 1,\*, Kamil Górecki <sup>2</sup> and Karolina Kaznowska-Opala <sup>3</sup>**


Received: 26 November 2018; Accepted: 4 January 2019; Published: 10 January 2019

**Abstract:** Water wells play an increasingly important role in providing water for the civilian population all over the world. Like other engineering structures, wells are subject to ageing processes resulting in degradation, which is observed as a reduction in hydraulic efficiency throughout their lifespan. To date, it has been found that the ageing process of a well is determined by a number of factors. The mathematical description of this process can be simplified. Drawing on Jacob's equation, this paper presents the course of the degradation process as a variable depending on operation time, well loss and flow rate. To apply the determined relationships in practice, simplifying assumptions were adopted, which make it possible to determine the moment of ageing compensations of the degradation processes. It was also demonstrated that the degradation process may be slowed down by the appropriate selection of initial operating parameters. The presented discussion highlights the significance of parameters *α*, *δ* and exponent *β*. The relation between hydraulic resistances in an aquifer and in the engineering structure is closely connected with these values. The presented arguments indicate that step drawdown tests provide the necessary information which allows tracking changes in the ageing processes occurring in the engineering structure. The analysis of the drawdown test results makes it possible to determine the moment when the necessary adjustments in the operating parameters of a water well should be performed. Eventually, it allows maintaining the high hydraulic efficiency of the intake and extending the lifespan of the well in accordance with the principle of sustainability.

**Keywords:** water well; hydraulic efficiency; degradation; engineering structure; well ageing; lifespan; well operation; water well management; sustainable efficiency

#### **1. Introduction**

Groundwater extraction is still growing all over the world. It is estimated that within half a century, i.e., between 1960 and 2010, groundwater use doubled [1]. The drilling of new groundwater wells is required. However, the number of boreholes is increasing not only due to population growth, but also due to a decrease in productivity of the existing assets. As well as other technical facilities wells are undergoing the process of degradation. As with any other technical facility, the current technical condition of a drilled well can be described with an efficiency index. Well efficiency is defined by Rorabaugh as the ratio of the theoretical drawdown computed by assuming that no turbulence is present to the drawdown in the well [2]. Walton defines the efficiency of a well as the ratio of the theoretical specific capacity to the actual specific capacity of the well [3]. Bierschenk concludes that well efficiency may be defined as the ratio of the theoretical drawdown to the measured drawdown inside the well. He also presents multiple examples of efficiency curves calculated using step-drawdown

tests carried out for wells located in the Middle East [4]. Therefore, well loss determines the hydraulic efficiency of a wellbore.

Well loss can be evaluated using the step-drawdown test. In Jacob's equation, well loss is described as the non-linear (non-laminar) flow regime component *CQ*<sup>2</sup> [5]. In some cases, changes of the drawdown can be explained by the Rorabaugh formula, where the power exponent *p* in the expression *CQ* is different from 2. Mackie, as cited by Atkinson et al. [6], reviewed the results of more than 20 carefully conducted step-drawdown tests of wells completed in fractured rock aquifers and concluded that most of the responses fell into one of the three categories where specific drawdown versus discharge rate is: linear, polynomial with power exponent equal to 2, or different from 2. The results of a study by Motyka and Wilk regarding the determination of the non-linear flow zone around several dozen wells drilled in fractured rocks indicate that the radius of turbulent flow zones are usually from 0.5 to 5.0 m, although in most cases they do not exceed1m[7].

Atkinson et al. concludes that high pumping rate moves the turbulent flow regime into the fissured aquifer, which is the reason for power exponent changes in Rorabaugh formula [8]. Klich et al. considered that the nature of the drawdown versus well-discharge curve depends on the range of well discharge [9]. Shekhar carried out the series of filed test for unconsolidated soils. He noticed that flow into the well screen is of a turbulent nature even for low discharge. The decrease in efficiency with increase of discharge can be considered from the perspective of increase in turbulence on account of the increase in discharge. As the drawdown increases the curvature of flow path increases leading to greater head loss. Therefore, well efficiency calculated can be regarded as a reflection of head loss on account of the laminar flow from the aquifer [10].

Much of the literature on the step-drawdown test focuses on the methods of determining aquifer resistance coefficient for linear water flow *B*, well resistance coefficient for non-linear flow *C* and *p* (if assumed not to be 2) parameters. Several methods are compiled by Kruseman and de Ridder [11]. One of the simplest graphical approaches is attributed to Jacob [5], Bruin and Hudson [12], Bierschenk and Wilson [13]. The parameters of the step-drawdown test can also be determined numerically. Miller and Webber [14] present an iterative method for solving the equation, while Labadie and Helweg developed a computer program for step-drawdown test data interpolation with the FASTEP procedure [15]. Avci proposes a method of analysis which calculates the aquifer and well loss from best fitting the log-log relationships of the difference of specific capacity and superposition of incremental pumping to the flow rate [16]. Wen et al. built the analytical model for all the types of curves considered for water well in confined aquifer. They noticed that all the curves at the wellbore approach the same asymptotic straight line in log s–log t scales [17]. Determining hydraulic characteristics of production wells from step-drawdown test data were also proposed by Jha et al. The characteristics of production wells such as aquifer loss coefficient, well loss coefficient, well specific capacity and well efficiency were determined by both the Genetic Algorithm optimization technique and the widely used graphical method. The developed computer programs also provide information about the condition of production wells and facilitate the construction of well characteristic curves [18].

Helweg proposes using General Well Function (GWF) to interpret the step-drawdown test. It allows more reliable results to be achieved for situations when the time that wells are continuously pumped varies greatly [19]. Kawecki presents a method for calculating total well loss based on the step drawdown test. In this method, the total loss can be estimated as a function of discharge rate. However, the real range radius of the borehole and also the specific storage of the confined aquifer are changeable at different stages of the step drawdown test [20]. An analytical solution for the constant pumping test in fissured porous media is presented by de Smedt. The solution is based on the dual-porosity approach with pseudo-steady state exchange between fissures and matrix. Proposed solution provides approach to analyze pumping test in fissured porous media characterized by double-porosity effects [21].

Singh proposes variable pumping replacing the step drawdown test, which does not require steady-state conditions at each step. This method allows the simultaneous estimation of both aquifer and well loss parameters [22]. The same author proposes a method for identifying head-loss from early drawdowns, where well loss could be evaluated from a short-duration pumping test in transient conditions. However, this method requires measurements in the pumped well and observation wells [23]. In contrast to this, Avci et al. proposed an analysis technique for interpreting transient step-drawdown test data. The method is based on taking the derivative of the drawdown with respect to time for the entire pumping test period, eliminating the constant well-loss terms. The conducted tests showed that the method allows the generation of the aquifer function *B*(*t*) for the entire duration of the step-drawdown test [24]. The choice of the model is the most difficult choice that the analyst of such a hydraulic test has to make since a wrong model can only lead to the wrong conclusions and failure of the borehole [25]. Although the physical model is well-known and widely used, there is a strong need to improve the techniques for estimating uncertainty associated with parameters derived from the interpretation of well tests [26].

The results of the step drawdown test could be applied to hydrogeology research. Dufresne presents the data used to develop a regional groundwater model which facilitates water planning and sustainability [27]. However, the step drawdown test is the most common tool to assess well performance and the hydraulic statement of the clogging. Many types of clogging processes are described in the book by Houben and Treskatis [28]. It involves a detailed discussion of the causes and effects of physical, chemical, electrochemical, and biochemical clogging. It also describes the methods and possibilities of rehabilitation of the well screen, gravel-pack and well-tube.

Van Beek presents the rehabilitation of mechanically clogged discharge wells. The research conducted so far has shown that the increase in flow velocity around a wellbore mobilizes fine particles, which are then transported through the porous structure. This may lead to the clogging of slots of the screen, thus reducing both porosity and permeability [29]. The mechanical clogging process in unconsolidated aquifers near the water supply wells is presented by de Zwart [30] and de Zwart et al. [31]. Blackwell et al. analyses particulate damage, which is often cited as a cause of permeability impairment around boreholes, resulting in a decline in well performance. Particulate movement and redistribution under borehole operating conditions have been assessed for a range of artificial and natural formations. Particulate damage cannot be effectively eliminated using normal development and rehabilitation techniques. This movement may occur at different stages in the life of the borehole, i.e., drilling, development and operational damage. Therefore, the operating regime has an important effect on well performance, which should be regularly monitored. If deterioration reaches approximately 10%, rehabilitation action can be initiated. Allowing deterioration to exceed 25% can greatly increase the cost of successful rehabilitation [32].

The example of chemical clogging including iron and iron-reducing bacteria is also presented by Hitchon [33]. The clogging of deep infiltration wells was also presented by de la Loma Gonzales, where several different types of chemical clogging and rehabilitation methods were compared and evaluated using a specific case of a well system [34]. The requirements for the economical construction and operation of drilled wells was also described by Treskatis et al [35]. These authors propose preventive maintenance for well rehabilitation as an essential guard against the development of non-soluble incrustations of Fe- and Mn-(hydr)oxides by back wash procedures or high-velocity, horizontal jetting techniques at an early stage of biochemical development aid to keep well efficiency at a high level. They also point out that surveillance of well efficiency and changes of drawdown can help to reduce the effort of chemicals for well rehabilitation. Van Beek concludes that well bore clogging may be prevented by regular intermittent abstraction. During rest, as there is no incoming flow, the particle accumulations will disintegrate. However, there is no clear picture of this disintegration process, an objective criterion for the rest time is lacking. It might be possible to shorten the rest period by reversing the flow direction with the abstraction flow velocity at least. Well bore clogging may also be counteracted by proper well development. Currently the criterion for a developed well is the absence of sand in the abstracted water. Actually, well development should be maximal until there is no further increase in specific capacity [36].

Houben et al. conclude that dramatic increases in head losses occur when clogging has reduced the effective porosity of the gravel pack by ~65%, the open area of the screen by ≥98%, and the casing diameter by ~50%. However, the clogging of the gravel pack is the most important object influencing well ageing. Moreover, obstructions in the gravel pack are much more difficult to remove. So, regular monitoring of well performance is needed, since processes in the gravel pack are difficult to track directly [37]. Therefore, low well efficiency contributes to the increased operating costs of wells. Helweg and Bengstone considered the economic approach the most profitable pumping rate in the whole life-span of the well, including of substitute well drilling [38]. Cost optimization throughout the life-span of a well is also the subject of a paper by Hurynowicz and Syczewa, who conclude that maintaining relatively low operating costs of a well depends on its appropriate design and operating conditions over several dozen years of its operation [39].

To sum up the above, well ageing causes well loss which is time-dependent. However, well loss is a naturally occurring phenomenon that depends on the hydraulic features of the aquifer, well construction, hydrogeochemistry, borehole drilling, well development and operation. The literature review presented above demonstrates that step-drawdown tests are typically used to estimate the hydrogeological parameters, well loss and well efficiency. Step-drawdown test results are also used to assess the current state and the effects of rehabilitation operations.

#### **2. Materials and Methods**

Water supply wells should be characterized by stable operation over a long period of time. They should therefore be operated with the highest possible hydraulic efficiency. After by van Tonder et al. [25], the rate at which a borehole can be pumped without lowering the water level below a set level is called "sustainable yield". Piscoppo and Summa discussed an experiment to estimate the pumping rate to ensure that the sustainable yield concept is realized [40]. At the beginning a step drawdown test was carried out, then the efficiency of the well was determined. Based on the results of the tests performed, the constant-head for the functioning of the well was defined. It was calculated that for a drawdown of around 12 m, the well efficiency is greater than 75%. The well was equipped with suitable submersible pump unit, then it was monitored on-line by one year. The pumping rate was defined as the discharge rate that will not cause the water level in the well to drop below a prescribed limit. In conclusions, authors emphasize that in some cases the constant-head pumping can be an alternative method of well management. Analysis of the recorded data changes shown that the discharge rate of the well trend with time was similar to that of the springs' hydrograms. Moreover, the water volume extracted did not exceed the recharge [40]. Subsequently, constant head model was compared with a constant flow rate model by Baiocchi et al. Based on numerical modelling, the authors concluded that constant flow rate model could be particularly useful when the problem is one of determining the sustainable yield of a single well from aquifers with low hydraulic diffusivity and when an extensive monitoring of the aquifer is not economically viable [41]. Unfortunately, the researchers did not evaluate whether both models had an impact on the well efficiency at the end of the research period.

In contrast to examples mentioned above, the authors of this paper have developed a concept of well operation while maintaining its constant efficiency, which could be also called "sustainable efficiency". This concept required the adoption of certain initial assumptions and was then determined theoretically using formulas evaluating the hydraulic efficiency of the well, taking into account the elements of the optimization calculus and also the ageing function of technical objects. Literature on hydrogeological issues states that the efficiency of a well can be described using an expression involving a rational function in the form proposed by Bierschenk [2]:

$$\begin{array}{l} \eta = \frac{s\_1}{s} = \frac{s\_1}{s\_1 + s\_2} = \frac{BQ}{\frac{BQ}{BQ} + CQ^2} = \frac{1}{1 + \frac{C}{B}Q} \\ \eta = \frac{1}{1 + aQ}, for \; a = \frac{C}{B} \end{array} \tag{1}$$

where:

*η* hydraulic efficiency [-];

*s*<sup>1</sup> aquifer loss [L];

*s* drawdown observed in the pumping well [L];

*s*<sup>2</sup> well loss in the well-screen adjacent zone [L];

*Q* volumetric flow rate [L3/T];

*B* aquifer resistance coefficient for laminar water flow [T/L2];

*C* well resistance coefficient for turbulent flow [T2/L5];

*α* parameter describing the prevalence of turbulent flow [T/L3].

This concept is presented in a schematic Figure 1.

**Figure 1.** Head losses in pumped well [4]; where: *Rt*—distance from center of well to effective point formation where transition from laminar to turbulent flow takes place.

At the beginning of the consideration presented below, it was assumed that the state of knowledge is sufficient to be used to select the optimum operating parameters of new wells and also to perform ongoing adjustments to counteract the intensive ageing of operated wells. So far, this subject matter has not been discussed in detail in the literature. For this reason, the authors present the results of their own research and considerations regarding the dynamics of the ageing process and the possibilities of reducing its rate. Below the article presents a process of degradation of a water well as an engineering structure. As previously mentioned, well efficiency function is time dependent. It means that a degradation process can be described as efficiency function in time. Therefore, well ageing is derivative from well efficiency in time.

For this purpose, the initial assumptions relating to the ageing process of the well could be defined by taking into account the state of knowledge presented in previous section:

**Assumption 1.** *Aquifer resistances are slow-variable. This assumption is valid for water supply wells where hydrogeological conditions do not change over time. In relation to drainage wells, the condition is usually not applicable due to the lowering of the water level and the relatively short life-span of the wells.*

**Assumption 2.** *The value of well loss (s2) should not significantly exceed the value of aquifer loss (s1). This assumption is usually relevant to new wells where the well-loss is insignificant in comparison to total drawdown. However, such criteria is also applicable to other wells in good hydraulic condition.*

Assumption 1 assumes that aquifer resistances are slow variables. This means that . *B*(*t*) ≈ 0 in comparison to the rates of changes in the well resistance parameters (*C*) or its discharge (Q)— . *<sup>B</sup>*(*t*) . *<sup>C</sup>*(*t*) (*or* . *Q*(*t*)) (accurate to the consistency factor of the dimension of a physical quantity). Therefore, if the said coefficients depend on time in a way that: *C* = *f*(*t*), *B* = *const*(*t*), *α* = *<sup>C</sup> <sup>B</sup>* = *g*(*t*) *and Q* = *h*(*t*), then well efficiency will also be a function of time *η* = *η*(*t*). The well efficiency function can be expanded into a Taylor series around moment *t* = *t*<sup>0</sup> as:

$$\eta(t) = \eta(t\_0) + \sum\_{n=1}^{k} \frac{1}{n!} \frac{d^{(n)}\eta}{dt^{(n)}} \bigg|\_{t=t\_0} (t - t\_0)^n \tag{2}$$

or by writing the above equation as follows:

$$\eta(t) = \eta(t\_0) + \frac{1}{1!} \dot{\eta}(t\_0) \cdot (t - t\_0) + \frac{1}{2!} \ddot{\eta}(t\_0) \cdot (t - t\_0)^2 + \frac{1}{3!} \dddot{\eta}(t\_0) \cdot (t - t\_0)^3 + \dots \tag{3}$$

Increases in the relevant variables may be used here—the increase in the independent variable (time) Δ*t* = *t* − *t*<sup>0</sup> and in the dependent variable (well efficiency) Δ*η* = *η* − *η*0. Then, the following relationship is present:

$$
\Delta \eta(t) = \dot{\eta}(t\_0) \cdot \Delta t + \frac{1}{2} \ddot{\eta}(t\_0) \cdot \Delta t^2 + \frac{1}{6} \dddot{\eta}(t\_0) \cdot \Delta t^3 + \dots \tag{4}
$$

The time derivative of well efficiency is:

$$\dot{\eta}(t) = \frac{d}{dt}(\frac{1}{1+aQ}) = \frac{-1}{(1+aQ)^2}(\dot{a}Q + a\dot{Q})\dot{\eta}(t) = -(\dot{a}Q + a\dot{Q})\cdot\eta^2\tag{5}$$

The second-order time derivative can be described by the expression:

$$\begin{aligned} \bar{\eta}(t) &= \frac{d}{dt}(\dot{\eta}) = (-1) \frac{(\dot{a}Q + a\dot{Q})^{\prime} \cdot (1 + aQ)^{2} - [(1 + Q)^{2}]^{\prime} \cdot (\dot{a}Q + a\dot{Q})}{[(1 + aQ)^{2}]^{2}} \\ \ddot{\eta}(t) &= (-1) \frac{(\ddot{a}Q + \dot{a}\dot{Q} + \dot{a}Q + a\dot{Q}) \cdot (1 + aQ)^{2} - 2(1 + aQ) \cdot (\dot{a}Q + a\dot{Q}) \cdot (\dot{a}Q + a\dot{Q})}{(1 + aQ)^{4}} \\ \ddot{\eta}(t) &= (-1) \frac{(\ddot{a}Q + 2\dot{a}Q + a\dot{Q}) \cdot (1 + aQ)^{2} - 2(1 + aQ) \cdot (\dot{a}Q + a\dot{Q})}{\left(1 + aQ\right)^{4}} \\ \ddot{\eta}(t) &= (-1) \frac{[(\ddot{a}Q + 2\dot{a}Q + a\dot{Q}) \cdot (1 + aQ) - 2(\dot{a}Q + a\dot{Q})^{2}] \cdot (1 + aQ)}{\left(1 + aQ\right)^{4}} \\ \ddot{\eta}(t) &= \frac{2(\ddot{a}Q + a\dot{Q})^{2} - (\ddot{a}Q + 2\dot{a}\dot{Q} + a\ddot{Q}) \cdot (1 + aQ)}{\left(1 + aQ\right)^{3}} \end{aligned} \tag{6}$$

The above can be expressed in a simplified manner as:

$$\begin{aligned} \bar{\eta}(t) &= \frac{2}{\eta} (\dot{\eta})^2 - \eta^2 (\ddot{a}Q + 2\dot{a}\dot{Q} + a\ddot{Q}) \\ \ddot{\eta}(t) &= \frac{2}{\eta} (\dot{\eta})^2 - \eta^2 \frac{d^2}{dt^2} (aQ) \end{aligned} \tag{7}$$

In Assumption 2 the value of well loss (*s*2) should not significantly exceed the value of aquifer loss (*s*1). This means that: *<sup>s</sup>*<sup>2</sup>

$$\begin{array}{c} \frac{\nu\_2}{s\_1} \ll 1\\ \frac{s\_2}{s\_1} = \frac{\mathbb{C}Q^2}{\mathbb{B}Q} = \frac{\mathbb{C}}{\mathbb{B}}Q = \mathbb{a}Q \ll 1 \end{array} \tag{8}$$

Therefore, the rate of change is small in comparison with other elements in the expression for the second-order time derivative of well efficiency:

$$\begin{array}{c} \frac{d^2}{dt^2}(\frac{s\_2}{s\_1}) = \frac{d^2}{dt^2}(aQ) \approx 0\\ \ddot{\eta}(t) = \frac{2}{\eta}(\dot{\eta})^2 - \eta^2 \frac{d^2}{dt^2}(aQ) \approx \frac{2}{\eta}(\dot{\eta})^2\\ \ddot{\eta}(t) \approx \frac{2}{\eta}(\dot{\eta})^2 \end{array} \tag{9}$$

Including in the expression the first-order derivative in the form: . *η*(*t*) = −( . *αQ* + *α* . *<sup>Q</sup>*)·*η*<sup>2</sup> provides the clear form: ..

$$\begin{aligned} \bar{\eta}(t) &= \frac{2}{\eta} \left[ (-1)(\dot{a}Q + a\dot{Q}) \cdot \eta^2 \right]^2 \\ \bar{\eta}(t) &= 2(\dot{a}Q + a\dot{Q})^2 \cdot \eta^3 \end{aligned} \tag{10}$$

The above considerations can be simplified by applying approximations in the solution, in which higher-order terms are negligible.

#### **3. Results**

The above equations describe the function of well degradation as a derivative of efficiency in time. However, it can be used to determinate the rate of ageing as well as to counteract this process. Issues can be considered in a simplified form as first-order approximation (kinetic), or in expansions of second-order approximations (dynamic). The higher-order terms, starting from the third order may be considered negligible. Changes in the function of the state of the well are presented in two independent variants, applying the principle of caeteris paribus. The first variant describes the ageing of the object and the decrease in hydraulic efficiency. The second variant describes the recovery of the lost efficiency of the well by changing the operational parameters of the well.

#### *3.1. First-Order Approximation*

In this case, higher-order terms are negligible, starting from the second order. In such a case the expression expanding the Taylor series for a function to a power series is identical to the expression for the differential of the function (adjacent to point (*t*0; *η*(*t*0)):

$$
\Delta \eta(t) \approx \dot{\eta}(t\_0) \cdot \Delta t \tag{11}
$$

which, after inserting the time derivative value for moment *t*<sup>0</sup> equals:

$$
\Delta \eta = -(\dot{a}Q + a\dot{Q}) \cdot \eta\_0^2 \cdot \Delta t \tag{12}
$$

To ensure clarity, subscripts are omitted for derivatives of functions *Q* and *α* at the moment of time *t*0. However, this fact should be taken into account when performing numerical calculations.

Time ranges can be matched in such a way that only one of the parameters, *α* or *Q*, will change distinctively over time, which would allow it to provide a significant contribution to the expression for changes in well efficiency throughout its lifespan. This will be discussed below on two extreme cases, which are viable on an engineering scale.

*Case 1–Well Ageing*

In this case it was assumed that the flow rate is invariable or slow-variable, while well resistances in relation to aquifer resistances increase considerably. This causes the deepening of the dynamic water table while maintaining a constant flow rate. This leads to a considerable reduction in the hydraulic efficiency of the well. Therefore, the turbulent flow of the reservoir fluids migrating in the near well zone acquires a high significance with time. This is connected with well ageing, which is a result of sediment accumulation in the functional part of the screen, clogging of the porous or fissure area of the medium directly adjacent to the screen and even the appearance of particulates inside the screen of the pumping well.

In other words: . *<sup>Q</sup>* <sup>≈</sup> <sup>0</sup> *and* . *α* > 0. Therefore:

$$
\Delta \eta\_{\text{ag\text{e}}} = - (\dot{a}Q + \dot{a}\dot{Q}) \cdot \eta\_0^2 \cdot \Delta t\_{\text{ag\text{e}}} \Delta \eta\_{\text{ag\text{e}}} = - \dot{a}Q \cdot \eta\_0^2 \cdot \Delta t\_{\text{ag\text{e}}} \tag{13}
$$

As can be easily determined, the inequality Δ*ηage <* 0 over time will always be satisfied. *Case 2—Compensation of Ageing*

In this case the flow rate decreases with time and the well resistances do not increase significantly. This leads to the compensation of ageing in the operated pumping well by reducing the flow rate. For . *<sup>Q</sup>* <sup>&</sup>lt; <sup>0</sup> *and* . *α* ≈ 0 the following can be specified:

$$
\Delta \eta\_{cmp} = -(\dot{a}Q + a\dot{Q}) \cdot \eta\_1^2 \cdot \Delta t\_{cmp} \\
\Delta \eta\_{cmp} = -a\dot{Q} \cdot \eta\_1^2 \cdot \Delta t\_{cmp} \tag{14}
$$

Obviously, in the above case the inequality Δ*ηcmp >* 0 over time will always be satisfied.

It should be pointed out here that a well with an initial efficiency of *η*<sup>0</sup> after the expiry of the period of ageing Δ*tage* reaches a state with the hydraulic efficiency of *η*1. This state represents the initial state for the later compensation of the well loss value. After the compensation time Δ*tcmp* the well reaches an efficiency of *η*. The number of iterations depends on the manner of transition between the extreme states of the well desired by the user.

#### *3.2. Second-Order Approximation*

In this case, it was assumed that higher-order terms, starting from the third order, are negligible. Therefore, the power series expansion of the well efficiency function is as follows:

$$
\Delta \eta(t) \approx \dot{\eta}(t\_0) \cdot \Delta t + \frac{1}{2} \ddot{\eta}(t\_0) \cdot \Delta t^2 \tag{15}
$$

while ensuring that the time derivatives at the moment of time *t=t*<sup>0</sup> are equal to:

$$\begin{aligned} \dot{\eta}(t\_0) &= -(\dot{a}Q + a\dot{Q}) \cdot \eta\_0^2\\ \ddot{\eta}(t) &\approx \frac{2}{\eta\_0} [\dot{\eta}(t\_0)]^2 = 2(\dot{a}Q + a\dot{Q})^2 \cdot \eta\_0^3 \end{aligned} \tag{16}$$

Therefore, the relationship for the well efficiency change will be equal to:

$$\begin{aligned} \Delta \eta(t) &\approx -(\dot{a}Q + a\dot{Q}) \cdot \eta\_0^2 \cdot \Delta t + \frac{1}{2} 2(\dot{a}Q + a\dot{Q})^2 \cdot \eta\_0^3 \cdot \Delta t^2 \\ \Delta \eta(t) &\approx -(\dot{a}Q + a\dot{Q}) \cdot \eta\_0^2 \cdot \Delta t + (\dot{a}Q + a\dot{Q})^2 \cdot \eta\_0^3 \cdot \Delta t^2 \\ \Delta \eta(t) &\approx (\dot{a}Q + a\dot{Q}) \cdot \eta\_0^2 \cdot \Delta t [(\dot{a}Q + a\dot{Q}) \cdot \eta\_0 \cdot \Delta t - 1] \end{aligned} \tag{17}$$

The first element of the difference provided in square brackets could be significantly smaller than one. This is possible because the time derivative of the product . *αQ* + *α* . *Q* = *<sup>d</sup> dt*(*αQ*) for a short time (where <sup>Δ</sup>*<sup>t</sup>* <sup>→</sup> 0) is low or slow-variable. Therefore, for *<sup>d</sup> dt*(*αQ*) ≈ 0 or *αQ* ≈ *const* and Δ*t* → 0 an accurate expression for well efficiency change can be obtained at first-order approximation.

#### *Case* 1*—Well Ageing*

Well ageing occurs with an invariable flow rate and an increase in the drawdown observed in the pumping well through a significant increase in well resistances and insignificant increase in aquifer resistances.

For . *<sup>Q</sup>* <sup>≈</sup> <sup>0</sup> *and* . *α* > 0 the simplified form of relationship (17) applies:

$$
\Delta \eta\_{\text{agc}} = -\dot{a}Q \cdot \eta\_0^2 \cdot \Delta t\_{\text{agc}} [1 - \dot{a}Q \cdot \eta\_0 \cdot \Delta t\_{\text{agc}}] \tag{18}
$$

The above equation will be satisfied for the loss of hydraulic efficiency of a well (Δ*ηage <* 0) if the passage of time is expressed in a limited range, given as the strict inequality Δ*tage* < ( . *αQ*·*η*0) −1 .

#### *Case* 2*—Well Ageing Compensation*

Flow rate decreases in an appropriately short time, while well resistances do not significantly increase over this time. For . *<sup>Q</sup>* <sup>&</sup>lt; <sup>0</sup> *and* . *α* ≈ 0 it will be:

$$
\Delta \eta\_{cmp} = -a \dot{Q} \cdot \eta\_1^2 \cdot \Delta t\_{cmp} [1 - a \dot{Q} \cdot \eta\_1 \cdot \Delta t\_{cmp}] \tag{19}
$$

As can be easily noticed, expression (19) describes inequality Δ*ηcmp >* 0 with time and decreasing flow rate in the well over time. In other words, it is possible to recover at least a part of the hydraulic efficiency lost. As before, after time Δ*tage*, well efficiency decreased from *η*<sup>0</sup> to *η*1, and it constitutes the initial value for the compensation of the well loss through a reduction in the flow rate of the borehole at time Δ*tcmp*. The number of iterations depends on the manner of transition between the extreme states of the well desired by the user.

#### **4. Discussion**

The archived results can be used in practice to compensate for well ageing. Taking into account that equations are difficult to use in practice, they can be expressed in a more practical way. Particular emphasis is put on the first order-approximation. For this purpose, characteristic indicators for full compensation of ageing will be presented below. In addition, expressions defining the necessary intervals to be taken in order to carry out corrective action are presented. The usage of those indicators for the selected example will be presented in the following section.

#### *4.1. First Order Approximation*

#### 4.1.1. Compensation of Well Ageing

It follows from case 1 that efficiency *η*<sup>1</sup> is:

$$
\eta\_1 = \eta\_0 - \dot{\kappa} Q \cdot \eta\_0^2 \cdot \Delta t\_{a\%} \tag{20}
$$

while obviously always Δ*ηage* = *η*<sup>1</sup> − *η*<sup>0</sup> < 0. Case 2 indicates that the final efficiency *η* value is:

$$
\eta = \eta\_1 - a \dot{Q} \cdot \eta\_1^2 \cdot \Delta t\_{cmp} \tag{21}
$$

while this time always Δ*ηcmp* = *η* − *η*<sup>1</sup> > 0.

The insertion of the expression for efficiency *η*<sup>1</sup> results in a change in efficiency as a result of compensation as: .

$$\begin{array}{c} \Delta \eta\_{cmp} = -aQ \cdot \eta\_1^2 \cdot \Delta t\_{cmp} \\ \Delta \eta\_{cmp} = -a\dot{Q} \cdot \left(\eta\_0 - \dot{a}Q \cdot \eta\_0^2 \cdot \Delta t\_{açc}\right)^2 \cdot \Delta t\_{cmp} \\ \Delta \eta\_{cmp} = -a\dot{Q} \cdot \eta\_0^2 \cdot \Delta t\_{cmp} \left(1 - \dot{a}Q \cdot \eta\_0 \cdot \Delta t\_{açc}\right)^2 \end{array} \tag{22}$$

#### 4.1.2. Full Compensation of Well Ageing

It is readily visible that the total change of well efficiency equals:

$$
\Delta \eta = \eta - \eta\_0 = (\eta - \eta\_1) + (\eta\_1 - \eta\_0) \\
\Delta \eta = \Delta \eta\_{cump} + \Delta \eta\_{\text{age}} \tag{23}
$$

the sign of the total change (difference) in efficiency is determined by the size of the individual elements.

Assuming that time moves towards moment *t*0, i.e., *t* → *t*0, so the time range becomes infinitesimal (Δ*ti* → 0) then the decrease in well efficiency as a result of ageing also becomes infinitesimal (Δ*ηage* → 0 or *η*<sup>1</sup> → *η*0). Furthermore, an increase in well efficiency as a result of compensating the well loss will also be infinitesimal (Δ*ηcmp* → 0 or *η* → *η*1). As for the rate of each of the above changes, it will be given using the "0/0" indeterminate form. However, it follows from the principle of transitivity of implication that *η* → *η*<sup>0</sup> also occurs, therefore in general the transition of the function to the limit Δ*η* → 0 where Δ*t* → 0 must occur. This means that the equality of function limits will be satisfied.

$$\lim\_{\Delta t \to 0} \Delta \eta = \lim\_{\Delta t \to 0} \Delta \eta\_{\text{ag}\varepsilon} = \lim\_{\Delta t \to 0} \Delta \eta\_{\text{cusp}} = 0 \tag{24}$$

Equations (23) and (24) provide knowledge about the possibility of the following relationship:

$$
\Delta \eta = \Delta \eta\_{\text{age}} + \Delta \eta\_{\text{comp}} = 0 \tag{25}
$$

It is always satisfied not only for initial moments of time, but also independently of time, which occurs with a full compensation of the decrease in well efficiency.

The insertion of expressions (13) and (14) into (25) results in:

$$
\dot{a}Q \cdot \eta\_0^2 \cdot \Delta t\_{a\xi c} + a \dot{Q} \cdot \eta\_1^2 \cdot \Delta t\_{cmp} = 0\tag{26}
$$

Applying Leibniz's rules to the differential calculus for function composition (the chain rule), it can be specified that if the proposition *d f dx* <sup>=</sup> *d f dt dx dt* <sup>⇒</sup> *d f dx dx dt* <sup>=</sup> *d f dt* is true, then the time from derivatives in expression (26) can be easily erased to obtain a first-order homogenous differential equation depending on the variable parameters *α* and *Q*:

$$\frac{d\boldsymbol{\alpha}}{d\boldsymbol{Q}} \cdot \eta\_0^2 \cdot \Delta t\_{\text{g}\text{g}^\text{e}} + \frac{\boldsymbol{\alpha}}{\boldsymbol{Q}} \cdot \eta\_1^2 \cdot \Delta t\_{\text{cmp}} = 0 \tag{27}$$

The separation of variables makes it possible to express the formula in the following form:

$$\begin{array}{c} -\frac{d\boldsymbol{a}}{\boldsymbol{a}} = \frac{d\boldsymbol{Q}}{\boldsymbol{Q}} \cdot \left(\frac{\boldsymbol{\eta}\_{1}}{\boldsymbol{\eta}\_{0}}\right)^{2} \left(\frac{\Lambda t\_{\rm cvp}}{\Lambda t\_{\rm gv}}\right) \\\ -\frac{d\boldsymbol{a}}{\boldsymbol{a}} = \beta \frac{d\boldsymbol{Q}}{\boldsymbol{Q}} \end{array} \tag{28}$$

The integration of relationship (28) with the constant of integration written as ln*γ<sup>0</sup> = const* and assuming for further simplicity constancy of coefficient *β* will result:

$$\begin{array}{lcl} -\int \frac{da}{a} = \int \beta \frac{dQ}{Q} = \beta \int \frac{dQ}{Q} \\ -\ln a = \beta \ln Q + const \\ \ln a^{-1} = \ln Q^{\delta} + \ln \gamma\_0 = \ln \left(\gamma\_0 Q^{\beta}\right) \end{array} \tag{29}$$

After introducing coefficient *δ* = *γ*−<sup>1</sup> <sup>0</sup> the above expression can be written in an antilogarithm (exponential) formula:

$$\begin{aligned} \alpha^{-1} &= \gamma\_0 Q^{\beta} \\ \gamma\_0^{-1} &= \delta = \alpha Q^{\beta} \end{aligned} \tag{30}$$

where exponent *β* of the power is equal to:

$$\beta = \left(\frac{\eta\_1}{\eta\_0}\right)^2 \left(\frac{\Delta t\_{cmp}}{\Delta t\_{ag\varepsilon}}\right) \tag{31}$$

Exponent *β* can assume different values, but some characteristic values can be listed. The following is a discussion of cases where: *β =* 0 and *β =* 1.

When *β =* 0, it should be expected that there is a strict inequality for the hydraulic efficiencies of the well at specific moments of time, i.e., *<sup>η</sup>*<sup>1</sup> <sup>&</sup>lt; *<sup>η</sup>*0. Therefore, their quotient is *<sup>η</sup>*<sup>1</sup> *<sup>η</sup>*<sup>0</sup> < 1 and the square of quotient is ( *<sup>η</sup>*<sup>1</sup> *η*0 ) <sup>2</sup> 1. Furthermore, it can be assumed that the time available for compensating well loss is shorter than the well ageing time. It is therefore realistic to assume that the relationship Δ*tcmp* < Δ*tage* will be true each time. This means that the exponent in Equation (31) will be equal to *β* << 1. In the context of mathematical analysis, it proceeds in its limit to the value of *β* → 0. Hence, expression (30) must assume the form:

$$\begin{aligned} \lim\_{\beta \to 0} \delta = \lim\_{\beta \to 0} (aQ^{\beta}) = \lim\_{\beta \to 0} (aQ^{0}) = \lim\_{\beta \to 0} a = a\\ \lim\_{\beta \to 0} \delta = a = \frac{\gamma}{B} = \delta\_0 \end{aligned} \tag{32}$$

Therefore, when the user of a water supply well is determined to achieve full compensation of the well loss in a wellbore used over a span of many years, then a certain stable state between the resistances of the aquifer and the resistances of the well should be established for the water flow in the aquifer. Indeed, as provided by Equation (32):

$$
\mathbb{C} = \mathfrak{a} \cdot \mathcal{B} \text{ or } \mathbb{C} = \delta\_0 \cdot \mathcal{B} \text{ when } \mathcal{B} \approx 0 \tag{33}
$$

The delta (*δ*0) parameter does depend on the property of the area directly adjacent to the functional part of the well screen and the permeability of the rock-soil medium conducting the fluid, as well as the type of the flowing reservoir medium.

Taking into account Equations (1) and (33), the hydraulic efficiency of the well can be equal to:

$$\eta = \frac{1}{1 + \alpha Q} = \frac{1}{1 + \delta\_0 Q} \text{ when } \beta \approx 0 \tag{34}$$

In other words, the ageing and compensation time ranges assumed for the well are similar in value (Δ*tage* ≈ Δ*tcmp*) and the well efficiencies compared to the limits of the structure's ageing period range are almost identical (*η*<sup>0</sup> ≈ *η*1), then the exponent from Equation (31) equals *β* ≈ 1. Obviously, this case will occur for initial moments—in situations where Δ*tage* → 0 and Δ*tcmp* → 0, i.e., where *η*<sup>1</sup> → *η*<sup>0</sup> (or Δ*ηage* → 0). Therefore, the limit:

$$\begin{array}{c} \lim\_{\beta \to 1} \delta = \lim\_{\beta \to 1} (aQ^{\beta}) = \lim\_{\beta \to 1} (aQ^{1}) = \lim\_{\beta \to 1} aQ = aQ \\ \lim\_{\beta \to 1} \delta = a \cdot Q = \frac{c}{B} \cdot Q = \delta\_{1} \end{array} \tag{35}$$

Then:

$$
\delta\_1 = \mathbf{a} \cdot \mathbf{Q} = \frac{\mathbf{C}}{B} \cdot \mathbf{Q} = \frac{\mathbf{C}Q^2}{BQ} \text{ or } \delta\_1 = \frac{s\_2}{s\_1} \text{ when } \not\gg \approx 1\tag{36}
$$

Parameter *δ*<sup>1</sup> for the initial moment under the conditions described above defines the mutual relationships between well loss and aquifer depression (or generally for a pool of any reservoir medium).

The hydraulic efficiency of a water supply well in accordance with the initial Jacob's formula, Equation (1), is then equal to:

$$\eta = \frac{s\_1}{s\_1 + s\_2} = \frac{s\_1}{s\_1 + \delta\_1 \cdot s\_1} = \frac{1}{1 + \delta\_1} \quad \text{when} \quad \beta \approx 1 \tag{37}$$

A table listing was prepared to summarize the discussion on the full compensation of well loss in the operated pumping well (Table 1).


**Table 1.** A listing of the calculation results for the full compensation of ageing-related well efficiency decrease, relationship of the type *<sup>δ</sup>* <sup>=</sup> *<sup>α</sup>*·*Q<sup>β</sup>* for <sup>Δ</sup>*<sup>η</sup>* <sup>=</sup> <sup>Δ</sup>*ηage* <sup>+</sup> <sup>Δ</sup>*ηcmp* <sup>=</sup> 0.

#### 4.1.3. Compensation Time

From Equations (30) and (31) the time required for the full compensation of hydraulic efficiency of the water supply well can be calculated. The following expression can be written:

$$\beta = \frac{\ln \frac{\delta}{a}}{\ln Q} = \left(\frac{\eta\_1}{\eta\_0}\right)^2 \left(\frac{\Delta t\_{cmp}}{\Delta t\_{q\xi^c}}\right) \tag{38}$$

Therefore, compensation time can be calculated using the formula:

$$
\Delta t\_{cmp} = \frac{\ln \frac{\delta}{a}}{\ln Q} \cdot \left(\frac{\eta\_0}{\eta\_1}\right)^2 \cdot \Delta t\_{ags} \tag{39}
$$

#### *4.2. Second Order Approximation*

Omitting the second element of the factor provided in square brackets in Equations (18) and (19) makes it possible to obtain a first-order approximation in such special cases. Writing down further equations would blur clarity due to complicated notation. At this stage Equations (18) and (19) provide information which allows taking action aimed at improving the hydraulic efficiency of the well.

The above discussion clearly demonstrates that periodical step drawdown tests play a fundamental role in the assessment of hydraulic state and also the potential actions aimed at extending the period of intake operation at a high level of hydraulic efficiency. However, the increase in well loss could be observed during regular or occasional inspections of operating parameters, i.e., drawdown and flow rate. In practice, compensation of the loss of hydraulic efficiency is achieved by methods such as valve throttling, switching the type series of pump units or, alternatively, regulating the rotational speed of the pump impeller.

Under conditions of a full compensation of wellbore ageing, parameter *δ*<sup>0</sup> is numerically related to the flow of the reservoir medium for which the flow rate of a water supply well *Q*<sup>0</sup> = (*C*/*B*) −1 results in the hydraulic efficiency of the well equal to 0.5 (50%), i.e., where well loss (*s*2) is equal to aquifer loss (*s*1). An increase in the value of parameter *δ*<sup>0</sup> means that the well reaches the same level of hydraulic efficiency with a lower volumetric flow rate. This is combined with a progressive increase in the resistance coefficient value of the well in turbulent flow at the well-screen near zone in relation to the resistance coefficient of laminar flow in the aquifer. In other words, parameter *δ*<sup>0</sup> may constitute a reliable indicator of well condition.

For a sufficiently long period of operation of the well the exponent is *β* = 0. Then for factor *δ* = *δ*<sup>0</sup> = *α* = *C*/*B* the time required to compensate well loss Δ*tcmp* = 0, i.e., in the conditions of well operation, is Δ*tcmp* << Δ*tage*. This means that for the long and intensive operation of water supply wells, flow rate reductions and compensation of hydraulic efficiency loss must not be delayed. For exponent *β* = 1, i.e., at the initial stage of the lifespan of a water supply well, factor *δ* = *δ*<sup>1</sup> = *s*2/*s*<sup>1</sup> = *αQ* ensures compensation time at a level equal to the time of ageing, i.e., operation of the structure—Δ*tcmp* ≈ Δ*tage*, provided that the efficiencies assume similar values (*η<sup>0</sup>* ≈ *η*1). In other words, the time available for taking corrective actions for the well becomes extended and is similar in its order of magnitude to the time of operation of the wellbore.

Furthermore, the higher the initial efficiency of the well (in relation to that expected after the period of operation of the well), the later it will be possible to begin compensating the loss of hydraulic efficiency by regulation through discharge, e.g., in an active or passive manner. On the other hand, excessively intensive pumping from the well will result in reducing the time available for compensating the effect of well loss. The well loss should increase with the clogging of the well-screen zone, manifested in the increase in the resistance coefficient value *C*, i.e., with the turbulent flow in the aquifer. Ultimately, this manifests itself in a decrease in the hydraulic efficiency curve of the well. Moreover, greater curvature of the graph is observed.

#### *4.3. Case Study*

Figures 2 and 3 present an example of using the solution presented above. Initially, Well 1 (Figure 2) operated at the flow rate of 12.5 × <sup>10</sup>−<sup>3</sup> m3/s, however parameter *<sup>α</sup>*−<sup>1</sup> <sup>=</sup> *B/C* <sup>=</sup> *<sup>Q</sup>*<sup>0</sup> = 37.4 × <sup>10</sup>−<sup>3</sup> m3/s. After 10 years of uninterrupted operation there was an unforeseen failure connected with the exposure of the pumping unit. A pumping test demonstrated a well loss of up to 80% of total drawdown. To protect the new pump unit, the pumping rate was limited to 5.6 × <sup>10</sup>−<sup>3</sup> m3/s, which corresponded to the hydraulic efficiency of the well of 30%. The conducted rehabilitation procedures failed to bring the expected results. The tests demonstrated that the final low efficiency of the well was caused by the relatively low initial well efficiency and the lack of efficiency compensation during operation. It has led to gravel-pack clogging and particulate damage. Finally, the value *Q*<sup>0</sup> decreased to 2.4 × <sup>10</sup>−<sup>3</sup> <sup>m</sup>3/s, which corresponds to an efficiency of 50%.

**Figure 2.** The hydraulic efficiency of the well without flow rate compensation: *t*0—at the beginning of operation time; *t*1—at the end of operation time.

**Figure 3.** The hydraulic efficiency of the well with flow rate compensation: *t*0—the working point at the beginning of operation; *t*1—the working point before first flow rate compensation; *t*1—the working point after first flow rate compensation; *t*2—the working point before second flow rate compensation; *t*2—the working point after second flow rate compensation; *η*0—the origin efficiency curve; *η*2—the current efficiency curve.

To reduce operating costs, the company conducting groundwater extraction constructed a replacement well. The methods presented in this article were used during its operation (Figure 3). Initially, the maintenance of an efficiency of 80% was established, which corresponded to the flow rate of 9.4 × <sup>10</sup>−<sup>3</sup> m3/s, however parameter *<sup>Q</sup>*<sup>0</sup> = 37.1 × <sup>10</sup>−<sup>3</sup> m3/s. After a further eight years the reduced flow rate was 5.8 × <sup>10</sup>−<sup>3</sup> <sup>m</sup>3/s with a constant efficiency. In this time, the value *<sup>Q</sup>*<sup>0</sup> decreased to 23.9 × <sup>10</sup>−<sup>3</sup> <sup>m</sup>3/s, which corresponds to an efficiency of 50%.

It should be noted that flow rate was significantly reduced for both wells. Currently, Well 1 is used occasionally due to high operating costs, while the costs of water extraction from Well 2 are approximately 65% lower. In addition, Well 2 operates at an efficiency which guarantees its reliable operation in the long term. The efficiency values, which are marked green and red in Figures 2 and 3, were confirmed by step drawdown tests. In practice, efficiency compensation in the second case was achieved by valve throttling at the beginning of operation and replacing the pump unit type series at appropriate time [16].

#### **5. Conclusions**

The solution presented in this paper involves ageing compensation and maintaining as high a hydraulic efficiency as possible. It concerns water supply wells in which well-loss has a minor share in the drawdown generated by water pumping. So, the presented methodology is valid to both new and old wells in good hydraulic condition, where hydrogeological conditions do not generally change over time. The fundamental meaning in presented method has parameter *α*. This parameter is a relation *C* to *B*, and may constitute a reliable well condition index which is relevant to 0.5 (50%) efficiency of the well.

Unlike "sustainable yield", concept of efficiency, compensation is not constant-rate or constant-head method. The "sustainable efficiency" concept requires regulation of both flow rate and drawdown. In relation to aquifer recharge fluctuations, operational parameters could be temporarily increased when water larger resources are available. However, in general, operating parameters of the well are reduced over time. This approach allows for the maintenance of a constant efficiency of the well, as long as the water output is sufficient to cover the water demand. Then, rehabilitation of the well is needed. Taking into consideration compensation time, it should be expressed that regulation time is dependent on current efficiency of the well. When the initial efficiency is relatively high, the compensation can be postponed or delayed. Alternatively, when the initial efficiency is low the time to reducing operational parameters is shortened.

The methodology was applied in a groundwater extraction company for new Well 2, drilled as a substitute of a Well 1, where efficiency was dramatically reduced over 10 years of constant water extraction. The periodical efficiency reduction of Well 2 during similar times of operation forced both the pumping rate and drawdown reduction. Currently, it has resulted in a relatively low cost of water extraction. Furthermore, the relation *B* to *C* still remains on a relatively high-level value. This indicates that the well-loss is quite small and the particulate damage in the near-well zone was minimized. The upcoming rehabilitation of Well 2 will probably be much more effective in this case. It can also be noticed that theoretical considerations presented in this paper were verified by step-drawdown tests at the beginning and also at the end of the research period. So, the method can be applied as a practical solution in well management to prevent inefficiencies.

**Author Contributions:** K.P. and K.K.-O. conceived of the presented idea. K.G. developed the theory and performed the computations. K.P. verified the analytical methods, performed the measurements and field test data calculations. K.K.-O. contributed to the interpretation of the results. All authors discussed the results and contributed to the final manuscript.

**Funding:** This research received no external funding. The APC was funded by AGH University of Science and Technology, Krakow 30-059, Poland.

**Acknowledgments:** The authors are grateful for the constructive comments made by anonymous reviewers which helped to improve the manuscript.

**Conflicts of Interest:** The authors declare no conflict of interest.

#### **References**


© 2019 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (http://creativecommons.org/licenses/by/4.0/).
