*Article* **Modeling of Artificial Groundwater Recharge by Wells: A Model Stratified Porous Medium**

**Carlos Fuentes 1, Carlos Chávez 2,\*, Antonio Quevedo 1, Josué Trejo-Alonso <sup>2</sup> and Sebastián Fuentes <sup>2</sup>**


Received: 21 September 2020; Accepted: 8 October 2020; Published: 13 October 2020

**Abstract:** In recent years, groundwater levels have been decreasing due to the demand in agricultural and industrial activities, as well as the population that has grown exponentially in cities. One method of controlling the progressive lowering of the water table is the artificial recharge of water through wells. With this practice, it is possible to control the amount of water that enters the aquifer through field measurements. However, the construction of these wells is costly in some areas, in addition to the fact that most models only simulate the well as if it were a homogeneous profile and the base equations are restricted. In this work, the amount of infiltrated water by a well is modeled using a stratified media of the porous media methodology. The results obtained can help decision-making by evaluating the cost benefit of the construction of wells to a certain location for the recharge of aquifers.

**Keywords:** mathematical modeling; infiltration well; differential equations; porous medium; fractal conductivity model

#### **1. Introduction**

Infiltration wells are used to contribute to the evacuation of rains in urban areas and also as a mechanism to recharge aquifers in regions where they present an unsustainable abatement [1–4]. Their construction must be analyzed from several angles: objectives of artificial recharge, available technological options, chemical quality of the water, social factors, place, quantity of water to contribute, among others [5–9].

In the literature, several numerical and analytical solutions can be found to model the flow of water in the porous medium, however, the models present restrictions to estimating the properties of soils, in addition to considering the stratum of the soil well profile as a homogeneous medium [1,10–14].

The artificial recharge capacity in a well is measured as the amount of water that infiltrates the soil during a specific period of time, and varies depending on the number of strata in the soil in which it was built. In this way, if you want to know the amount of water that the entire well contributes, you must evaluate the infiltration rate in all the strata to have a better knowledge about the contributions to the aquifer and the behavior of the system as a whole.

The phenomenon of infiltration in porous media can be studied from the general principles of the conservation of mass and momentum. The equation that results from the application of the first principle is:

$$\frac{\partial \Theta}{\partial t} = -\nabla \cdot \overrightarrow{\mathbf{q}}.\tag{1}$$

Darcy's law generalized to partially saturated porous media is used as a dynamic equation [15]:

$$
\overrightarrow{\mathbf{q}} = -\mathsf{K}\nabla\mathsf{H},\tag{2}
$$

where H is the hydraulic potential and is the sum of the pressure potential (ψ) and the gravitational potential assimilated to the vertical coordinate (z) oriented, in this case, as positive upwards. The pressure potential is positive in the saturated zone and negative in the unsaturated zone, since it is agreed that the zero pressure corresponds to the atmospheric pressure; θ = θ (ψ) is the volumetric water content, also called moisture content, and is a function of the water pressure, θ (ψ) is known as the retention curve or soil moisture characteristic; <sup>→</sup> q = (qx, qy, qz) is the flow of water per unit of soil surface or Darcy flow, with its components in a rectangular system; (x, y, z) are the spatial coordinates in a rectangular or Cartesian system, t is time; ∇ is the gradient operator; K = K (ψ) is the hydraulic conductivity as a function of the water pressure.

Thus, the general equation of flow in a porous medium results from the combination of Equations (1) and (2):

$$\frac{\partial \Theta}{\partial \mathbf{t}} = \nabla \cdot [\mathbf{K}(\boldsymbol{\upmu}) \nabla(\boldsymbol{\upmu} + \mathbf{z})].\tag{3}$$

This equation presents two independent variables, θ and ψ, but since there is a relationship between them, the specific capacity defined as the slope of the retention curve is introduced. The chain rule is applied and the equation with the dependent variable pressure is established, known as the Richards equation [16]:

$$\mathbf{C}(\boldsymbol{\psi})\frac{\partial\boldsymbol{\psi}}{\partial\mathbf{t}} = \nabla \cdot \left[\mathbf{K}(\boldsymbol{\psi})\nabla\boldsymbol{\psi}\right] + \frac{\partial\mathbf{K}}{\partial\boldsymbol{\psi}}\frac{\partial\boldsymbol{\psi}}{\partial\mathbf{z}};\,\mathbf{C}(\boldsymbol{\psi}) = \frac{\partial\boldsymbol{\theta}}{\partial\boldsymbol{\psi}}.\tag{4}$$

In this work a methodology is presented to obtain the water infiltration rate by partially or totally filled artificial recharge wells. The equations have been adapted to be used in a homogeneous stratified medium, taking into account the soil characteristics of each strata in the profile of the well.

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

#### *2.1. The Richards and Kirchho*ff *Equations in Spherical and Cylindrical Coordinates*

In some problems the analysis is simplified if Equation (4) is written in cylindrical or spherical coordinates. The Richards equation [16] in cylindrical coordinates (r, ϕ, z) is as follows:

$$\mathbf{C}(\boldsymbol{\Psi})\frac{\partial\boldsymbol{\Psi}}{\partial\mathbf{t}} = \frac{1}{\mathbf{r}}\frac{\partial}{\partial\mathbf{r}}\Big[\mathbf{r}\mathbf{K}(\boldsymbol{\Psi})\frac{\partial\boldsymbol{\Psi}}{\partial\mathbf{r}}\Big] + \frac{1}{\mathbf{r}^2}\frac{\partial}{\partial\boldsymbol{\varrho}}\Big[\mathbf{K}(\boldsymbol{\Psi})\frac{\partial\boldsymbol{\Psi}}{\partial\boldsymbol{\varrho}}\Big] + \frac{\partial}{\partial\mathbf{z}}\Big[\mathbf{K}(\boldsymbol{\Psi})\frac{\partial\boldsymbol{\Psi}}{\partial\mathbf{z}}\Big] + \frac{\partial\mathbf{K}}{\partial\boldsymbol{\Psi}}\frac{\partial\boldsymbol{\Psi}}{\partial\mathbf{z}}\Big\} \tag{5}$$

where r is the radius and ϕ is the azimuth: r<sup>2</sup> = x<sup>2</sup> + y2, x = r cos ϕ, y = r sin ϕ.

In spherical coordinates (, ϑ, ϕ) the Richards equation is written as:

$$\begin{array}{l} \mathbf{C}(\boldsymbol{\Psi}) \frac{\partial \boldsymbol{\Psi}}{\partial \boldsymbol{\Psi}} = \frac{1}{\varrho^{2}} \frac{\partial}{\partial \boldsymbol{\Psi}} \Big[ \varrho^{2} \mathbf{K}(\boldsymbol{\Psi}) \frac{\partial \boldsymbol{\Psi}}{\partial \boldsymbol{\varrho}} \Big] + \frac{1}{\varrho^{2} \sin \boldsymbol{\Psi}} \frac{\partial}{\partial \boldsymbol{\Psi}} \Big[ \sin \boldsymbol{\vartheta} \mathbf{K}(\boldsymbol{\Psi}) \frac{\partial \boldsymbol{\Psi}}{\partial \boldsymbol{\Psi}} \Big] \\ + \frac{1}{\varrho^{2} \sin^{2} \boldsymbol{\Psi}} \frac{\partial}{\partial \boldsymbol{\Psi}} \Big[ \mathbf{K}(\boldsymbol{\Psi}) \frac{\partial \boldsymbol{\Psi}}{\partial \boldsymbol{\Psi}} \Big] + \frac{\partial \mathbf{K}}{\partial \boldsymbol{\Psi}} \frac{\partial \boldsymbol{\Psi}}{\partial \mathbf{x}} \end{array} \tag{6}$$

where is the radio, ϑ is the polar angle and ϕ is the azimuth: <sup>2</sup> = x2 + y2 + z2, x = sin ϑ cos ϕ, y = sin ϑ sin ϕ, z = cos ϑ.

In a symmetric well with respect to the z axis, the radius r takes this axis as its origin, and Equation (5) is very useful for the infiltration analysis when it is assumed that the pressure does not depend on the azimuth, that is, when the heterogeneity is presented by layers. In this case the equation simplifies to the following:

*Mathematics* **2020**, *8*, 1764

$$\mathbf{C}(\boldsymbol{\Psi})\frac{\partial\boldsymbol{\Psi}}{\partial\mathbf{t}} = \frac{1}{\mathbf{r}}\frac{\partial}{\partial\mathbf{r}}\bigg[\mathbf{r}\mathbf{K}(\boldsymbol{\Psi})\frac{\partial\boldsymbol{\Psi}}{\partial\mathbf{r}}\bigg] + \frac{\partial}{\partial\mathbf{z}}\bigg[\mathbf{K}(\boldsymbol{\Psi})\frac{\partial\boldsymbol{\Psi}}{\partial\mathbf{z}}\bigg] + \frac{\partial\mathbf{K}}{\partial\boldsymbol{\Psi}}\frac{\partial\boldsymbol{\Psi}}{\partial\mathbf{z}}\tag{7}$$

which has only two spatial coordinates (r, z).

The equation in spherical coordinates presents a very particular importance when, in the analysis of a problem, it is considered that the medium is homogeneous and isotropic, that is, when the phenomenon does not depend on either the colatitude or the azimuth:

$$\mathcal{K}(\boldsymbol{\psi})\frac{\partial\boldsymbol{\psi}}{\partial\mathbf{t}} = \frac{1}{\varrho^2}\frac{\partial}{\partial\varrho}\bigg|\varrho^2\mathcal{K}(\boldsymbol{\psi})\frac{\partial\boldsymbol{\psi}}{\partial\varrho}\bigg| + \frac{\partial\mathcal{K}}{\partial\boldsymbol{\psi}}\frac{\partial\boldsymbol{\psi}}{\partial\mathbf{z}}\tag{8}$$

in which only two spatial coordinates are presented (, z).

Furthermore, in some particular problems involving homogeneous porous media, the analysis is simplified if the potential Kirchhoff flow is defined by:

$$\Phi = \int\_{-\infty}^{\psi} \mathbf{K}(\overline{\Psi}) \mathrm{d}\overline{\Psi} = \int\_{\Theta\_{\mathrm{r}}}^{\Theta} \mathrm{D}(\overline{\Theta}) \mathrm{d}\overline{\Theta}, \tag{9}$$

from which it follows that:

is:

$$\frac{d\Phi}{d\psi} = \mathcal{K}(\psi); \frac{d\Phi}{d\theta} = \mathcal{D}(\theta), \tag{10}$$

where D (θ) is the hydraulic diffusivity, in analogy with the diffusion of gases, which is expressed as D(θ) = K(θ)/C(θ), considering, now, that both the hydraulic conductivity and the specific capacity are functions of the volumetric content moisture.

The water transfer equation in porous media as a dependent variable for moisture content, Equation (3), is as follows:

$$\frac{\partial \Theta}{\partial t} = \nabla \cdot \left[ \mathbf{D}(\theta) \nabla \theta \right] + \frac{\mathbf{d} \mathbf{K}}{\mathbf{d}\theta} \frac{\partial \Theta}{\partial \mathbf{z}'} \tag{11}$$

which presents the structure of a nonlinear Fokker–Planck equation [17], the linear version of which is widely known in diffusion problems.

In terms of the potential Kirchhoff flow, Equation (11) becomes:

$$\frac{1}{\mathrm{D}(\Phi)}\frac{\partial\Phi}{\partial\mathbf{t}}=\nabla^{2}\Phi+\frac{\mathrm{d}\mathbf{K}}{\mathrm{d}\Phi}\frac{\partial\Phi}{\partial\mathbf{z}}.\tag{12}$$

Kirchhoff's equation in cylindrical coordinates (r, ϕ, z), where r is the radius and ϕ is the azimuth,

$$\frac{1}{\mathrm{D}(\Phi)}\frac{\partial\Phi}{\partial\mathbf{t}} = \frac{1}{\mathrm{r}}\frac{\partial}{\partial\mathbf{r}}\Big(\mathrm{r}\frac{\partial\Phi}{\partial\mathbf{r}}\Big) + \frac{1}{\mathrm{r}^2}\frac{\partial^2\Phi}{\partial\varphi^2} + \frac{\partial^2\Phi}{\partial\mathbf{z}^2} + \frac{\mathrm{dK}}{\mathrm{d}\Phi}\frac{\partial\Phi}{\partial\mathbf{z}}.\tag{13}$$

In spherical coordinates (, ϑ, ϕ) the Kirchhoff equation is written as follows:

$$\begin{array}{ll} \frac{1}{\mathrm{D}(\Phi)} \frac{\partial \Phi}{\partial t} = \frac{1}{\varrho^2} \frac{\partial}{\partial \varrho} \Big( \varrho^2 \frac{\partial \Phi}{\partial \varrho} \Big) + \frac{1}{\varrho^2 \sin \vartheta} \frac{\partial}{\partial \vartheta} \Big( \sin \vartheta \frac{\partial \Phi}{\partial \vartheta} \Big) \\ + \frac{1}{\varrho^2 \sin^2 \vartheta} \frac{\partial^2 \Phi}{\partial \varphi^2} + \frac{d \mathrm{K}}{d \Phi} \frac{\partial \Phi}{\partial \mathbf{z}} \end{array} . \tag{14}$$

These formulations are only applicable in homogeneous media, and, if these are isotropic, they are written respectively in cylindrical coordinates as follows:

$$\frac{1}{\mathrm{D}(\Phi)}\frac{\partial\Phi}{\partial\mathbf{t}} = \frac{1}{\mathrm{r}}\frac{\partial}{\partial\mathbf{r}}\Big(\mathbf{r}\frac{\partial\Phi}{\partial\mathbf{r}}\Big) + \frac{\partial^2\Phi}{\partial\mathbf{z}^2} + \frac{\mathrm{d}\mathbf{K}}{\mathrm{d}\Phi}\frac{\partial\Phi}{\partial\mathbf{z}}\Big) \tag{15}$$

and in spherical coordinates:

$$\frac{1}{\mathrm{D}(\Phi)}\frac{\partial\Phi}{\partial t} = \frac{1}{\varrho^2}\frac{\partial}{\partial\varrho}\bigg(\varrho^2\frac{\partial\Phi}{\partial\varrho}\bigg) + \frac{\mathrm{d}K}{\mathrm{d}\Phi}\frac{\partial\Phi}{\partial\mathbf{z}}.\tag{16}$$

#### *2.2. The Hydrodynamic Characteristics of Porous Media*

In order to solve the mass or energy transfer equations of water in porous media, aside from specifying the limit conditions, it is necessary to know the hydrodynamic characteristics formed by the water retention curve θ (ψ) and the hydraulic conductivity curve either as a function of the water pressure, K(ψ), or as a function of the moisture content K(θ). The analysis is greatly simplified if these curves are represented with analytical functions.

The retention curve can be represented with the equation of van Genuchten [18]:

$$\Theta(\psi) = \left[1 + \left(\frac{\psi}{\psi\_{\rm d}}\right)^{\rm n}\right]^{-m},\tag{17}$$

where m > 0 and n > 0 are two shape parameters (dimensionless), ψ<sup>d</sup> is a characteristic value of the water pressure and Θ is the effective degree of saturation defined by:

$$
\Theta = \frac{\theta - \theta\_{\rm r}}{\theta\_{\rm s} - \theta\_{\rm r}},
\tag{18}
$$

in which θ<sup>r</sup> is the residual moisture content defined such that K(θr) = 0 and θ(ψ→-∞) = θ<sup>r</sup> [19]; θ<sup>s</sup> is the moisture content at saturation, assimilated to the total porosity of the soil (φ), when under saturation conditions no air is trapped in the interstices of the porous medium: θ<sup>s</sup> = φ. In general, θ<sup>r</sup> = 0 can be assumed [20].

A closed way to represent the conductivity curve can be obtained using prediction models of the same from the retention curve. In the literature we can find various works, but given the condition of the phenomenon we are studying, this work uses one of the fractal models proposed, calibrated and validated by Fuentes et al. [20]. The results found by [20] shows a better adjustment between observed and estimated data by the function given by:

$$\mathcal{K}(\Theta) = \mathcal{K}\_s \left[ \int\_0^{\Theta} \frac{\vartheta^{s-1} \mathrm{d}\vartheta}{\left| \psi(\vartheta) \right|^{2s}} / \int\_0^1 \frac{\vartheta^{s-1} \mathrm{d}\vartheta}{\left| \psi(\vartheta) \right|^{2s}} \right]^2,\tag{19}$$

where Ks is the hydraulic saturation conductivity and s = D/E, with D the fractal dimension of the porous medium and E = 3 the Euclid dimension of the physical space where the medium is embedded, related to the total porosity through the relation:

$$(1 - \phi)^{s} + \phi^{2s} = 1.\tag{20}$$

The introduction of Equation (17) in Equation (19) leads to the following equation to represent the hydraulic conductivity curve, accepting the relationship between the parameters as indicated:

$$\mathcal{K}(\Theta) = \mathcal{K}\_{\mathbf{s}} \Big[ 1 - \left( 1 - \Theta^{1/\mathbf{m}} \right)^{\mathrm{sm}} \Big]; 0 < \mathrm{sm} = 1 - 2\mathrm{s}/\mathrm{n} < 1. \tag{21}$$

The solution of the transfer equation in its different forms is generally numerical [21,22]. However, in some simplified cases, characteristics of the solution can be obtained analytically [23–25].

#### **3. Results and Discussion**

#### *3.1. Conceptual Model*

To analyze the infiltration in steady-state wells, it is necessary to write the Darcy flows in the radial and vertical directions:

$$\stackrel{\rightharpoonup}{\mathbf{q}}\_{\mathbf{r}} = -\mathbf{K}\_{\mathbf{s}} \frac{\partial \psi}{\partial \mathbf{r}} \mathbf{\hat{r}}\_{\mathbf{r}} \tag{22}$$

$$\stackrel{\rightharpoonup}{\mathbf{q}}\_x = - \Big( \mathbb{K}\_\mathbf{\hat{s}} \frac{\partial \Psi}{\partial \mathbf{z}} + \mathbb{K}\_\mathbf{\hat{s}} \Big) \hat{\mathbf{k}}\_\mathbf{\hat{s}} \Big. \tag{23}$$

where r and ˆ k are unitary vectors in the r and z directions respectively. ˆ

Flow through the wall and bottom of the well is defined by:

$$\mathbf{Q}\_{\mathbf{s}} = \int\_{\mathbf{A}\_{\mathbf{P}}} \overrightarrow{\mathbf{q}}\_{\mathbf{r}} \cdot \mathbf{dA}\_{\mathbf{P}} + \int\_{\mathbf{A}\_{\mathbf{b}}} \overrightarrow{\mathbf{q}}\_{\mathbf{z}} \cdot \mathbf{dA}\_{\mathbf{b}\prime} \tag{24}$$

where dAp y dAb are, respectively, the differential areas in the wall and at the bottom of the well defined by:

$$\mathbf{dA}\_{\rm P} = (2\pi \mathbf{R} \mathbf{dz})\mathbf{\hat{r}},\tag{25}$$

$$\mathbf{dA}\_{\mathbf{b}} = (2\pi \mathbf{r} \mathbf{D} \mathbf{r}) \mathbf{(-\hat{k})}.\tag{26}$$

Equation (24), considering Equations (22), (23), (25) and (26), is written as follows:

$$\mathbf{Q}\_{\mathbf{s}} = -2\pi \mathbf{R} \mathbf{K}\_{\mathbf{s}} \int\_{0}^{\mathbf{H}} \left. \frac{\partial \boldsymbol{\psi}}{\partial \mathbf{r}} \right|\_{\mathbf{r}=\mathbf{r}} \, \mathrm{d}\mathbf{z} + 2\pi \mathbf{K}\_{\mathbf{s}} \int\_{0}^{\mathbf{R}} \left. \frac{\partial \boldsymbol{\psi}}{\partial \mathbf{z}} \right|\_{\mathbf{z}=0} \, \mathrm{r} \mathrm{d}\mathbf{r} + \pi \mathbf{K}\_{\mathbf{s}} \mathrm{R}^{2} . \tag{27}$$

Introducing the dimensionless variables:

$$\mathbf{z}^\* = \frac{\mathbf{z}}{\mathbf{H}}; \mathbf{r}^\* = \frac{\mathbf{r}}{\mathbf{R}}; \psi^\* = \frac{\psi}{\mathbf{H}}; \tag{28}$$

Equation (27) is written as follows:

$$\mathbf{Q}\_{\text{s}} = \mathbf{Q}\_{\text{o}} + \pi \mathbf{K}\_{\text{s}} \mathbf{R}^2; \; \mathbf{Q}\_{\text{o}} = \frac{2\pi \mathbf{K}\_{\text{s}} \mathbf{H}^2}{\mathbf{C}};\tag{29}$$

where C is a coefficient defined as:

$$\frac{1}{\mathbf{C}} = -\int\_0^1 \frac{\partial \boldsymbol{\upmu}^\*}{\partial \mathbf{r}^\*}\Big|\_{\mathbf{r}^\* = 1} \, \mathrm{d}\mathbf{z}^\* + \left(\frac{\mathrm{R}}{\mathrm{H}}\right)^2 \int\_0^1 \frac{\partial \boldsymbol{\upmu}^\*}{\partial \mathbf{z}^\*}\Big|\_{\mathbf{r}^\* = 0} \, \mathrm{r}^\* \mathrm{d}\mathbf{r}^\*.\tag{30}$$

To find this coefficient it is necessary to know ψ (r, z).

#### *3.2. The Glover Model*

According to Glover [26], in a first approximation, the steady-state pressure flow through an infiltration well in a homogeneous and isotropic porous medium can be described with Laplace's equation in spherical coordinates that describes the pressure in the absence of gravitational gradients. From Equation (8) we have:

$$\nabla^2 \Psi = \frac{1}{\varrho^2} \frac{\partial}{\partial \varrho} \bigg( \varrho^2 \frac{\partial \Psi}{\partial \varrho} \bigg) = 0,\tag{31}$$

which must be subject to border conditions:

$$
\Psi = \Psi\_{\mathbb{R}} ; \; \varrho = \mathbb{R} \tag{32}
$$

$$
\Psi = 0; \ \varrho \to \infty. \tag{33}
$$

Integration of Equation (31) leads to <sup>ψ</sup> = <sup>−</sup>c1−<sup>1</sup> + c2, where c1 and c2 are integration constants; the Equation (33) implies c2 = 0 and the Equation (32) c1 = −ψ<sup>R</sup> R, ergo ψ = ψ<sup>R</sup> (R/). The Darcy flux is <sup>q</sup> <sup>=</sup> <sup>−</sup> Ks <sup>∂</sup>ψ/∂ <sup>=</sup> Ks <sup>ψ</sup><sup>R</sup> (R/2), when <sup>=</sup> R, qR <sup>=</sup> Ks <sup>ψ</sup>R/R; the flow through the surface of the sphere of radius R is qo = 4πR2qR = 4πKsRψR; this flow from the point source in the center of the sphere is the variable of interest. Since ψRR = qo/4πKs, it is better to set the pressure variation around the source flow to continue with the Glover approach:

$$
\Psi = \frac{\mathbf{q}\_{\text{o}}}{4\pi \mathbf{K}\_{\text{s}}\varrho}.\tag{34}
$$

If h represents the position of the center of the sphere from the base, then the spherical coordinate () and the cylindrical coordinate (r) are related by:

$$
\varrho = \sqrt{\mathbf{r}^2 + \left(\mathbf{z} - \mathbf{h}\right)^2}.\tag{35}
$$

The pressure in terms of the cylindrical coordinates is obtained by introducing Equation (35) into Equation (34):

$$\Psi = \frac{\mathbf{q}\_{\text{o}}}{4\pi\mathbf{K}\_{\text{s}}\sqrt{\mathbf{r}^{2} + \left(\mathbf{z} - \mathbf{h}\right)^{2}}} \cdot \tag{36}$$

To provide a series of point sources whose magnitude increases with depth, an expression similar to that originally proposed by Glover, we have:

$$\mathbf{dq}\_{\rm o} = \mathbf{B}(\mathbf{h}\_{\rm \tilde{c}} - \mathbf{h})\mathbf{dh}\_{\rm \tilde{s}} \tag{37}$$

where B is a parameter to be determined and hc defines the range of the sources ho ≤ h ≤ hc and sinks hc < h ≤ hs.

The total flow is found by integrating Equation (37):

$$\mathbf{Q}\_{\rm o} = \underset{\mathbf{H}\_{\rm o}}{\text{B}} \stackrel{\text{h}\_{\rm o}}{\text{C}} (\mathbf{h}\_{\rm c} - \mathbf{h}) \text{dh} = \frac{1}{2} \text{BH}^2 \Big[ \left( \mathbf{h}\_{\rm c}^\* - \mathbf{h}\_{\rm o}^\* \right)^2 - \left( \mathbf{h}\_{\rm c}^\* - \mathbf{h}\_{\rm s}^\* \right)^2 \Big] \tag{38}$$

hence parameter B is deduced:

$$\mathbf{B} = \frac{2\mathbf{Q}\_0}{\mathbf{H}^2 \left[ \left( \mathbf{h}\_\mathrm{c}^\* - \mathbf{h}\_\mathrm{o}^\* \right)^2 - \left( \mathbf{h}\_\mathrm{c}^\* - \mathbf{h}\_\mathrm{s}^\* \right)^2 \right]} \tag{39}$$

where h\* = h/H for all subscripts.

From Equations (36), (37) and (39) we have:

$$\mathrm{d}\Psi = \frac{\mathrm{Q}\_{\mathrm{o}} \mathrm{(h\_{c} - h)} }{2\pi \mathrm{K}\_{\mathrm{e}} \mathrm{H}^{2} \left[ \left( \mathrm{h\_{c}^{\*} - \mathrm{h\_{o}} \right)^{2} - \left( \mathrm{h\_{c}^{\*} - \mathrm{h\_{s}^{\*}} \right)^{2} \right]} \sqrt{\mathrm{r}^{2} + \left( \mathrm{z} - \mathrm{h} \right)^{2}}} \mathrm{d}\mathrm{h} \tag{40}$$

the integration of which leads to:

$$\Psi = \frac{\mathbf{Q\_o}}{2\pi \mathbf{K\_s} \mathbf{H^2} \left[ \left( \mathbf{h\_c^\* - h\_o^\*} \right)^2 - \left( \mathbf{h\_c^\* - h\_s^\*} \right)^2 \right]} \left( \mathbf{h\_c} - \mathbf{z} \right) \mathbf{a} \sinh \left( \frac{\mathbf{z} - \mathbf{h}}{\mathbf{r}} \right) + \sqrt{\mathbf{r}^2 + \left( \mathbf{z} - \mathbf{H} \right)^2} \bigg|\_{\mathbf{h} = \mathbf{h\_s}}^{\mathbf{h} - \mathbf{h\_o}}, \tag{41}$$

ergo:

$$\Psi = \frac{\mathbf{Q}\_o}{2\pi \mathbf{K}\_s \mathbf{H}^2} \frac{\left[ \frac{\left( \mathbf{h}\_c - \mathbf{z} \right) \mathbf{as} \text{sinh} \left( \frac{\mathbf{z} - \mathbf{h}\_o}{\mathbf{r}} \right) - \left( \mathbf{h}\_c - \mathbf{z} \right) \mathbf{as} \text{sinh} \left( \frac{\mathbf{z} - \mathbf{h}\_o}{\mathbf{r}} \right)}{\sqrt{\mathbf{r}^2 + \left( \mathbf{z} - \mathbf{h}\_o \right)^2} - \sqrt{\mathbf{r}^2 + \left( \mathbf{z} - \mathbf{h}\_s \right)^2}} \right]}{\left[ \left( \mathbf{h}\_c^\* - \mathbf{h}\_o^\* \right)^2 - \left( \mathbf{h}\_c^\* - \mathbf{h}\_s^\* \right)^2 \right]}. \tag{42}$$

At the point on the boundary (r, z) = (R, 0) we have ψ = H, which allows obtaining the expression of the flow, Equation (29):

$$\mathbf{Q}\_{\rm o} = \frac{2\pi \mathbf{K}\_{\rm s} \mathbf{H}^2}{\mathbf{C}},\tag{43}$$

where the form coefficient is defined by:

$$\mathbf{C} = \frac{\mathbf{h}\_{\mathrm{c}}^{\*} \Big[ \mathrm{asinh} \Big( \frac{\mathbf{h}\_{\mathrm{c}}^{\*} \mathbf{h}\_{\mathrm{s}}^{\*} \Big) - \mathrm{asimh} \Big( \frac{\mathbf{h}\_{\mathrm{c}}^{\*} \mathbf{h}\_{\mathrm{o}}^{\*} \Big) \Big] + \sqrt{\left( \frac{\mathbf{R}}{\mathrm{H}} \right)^{2} + \mathrm{h}\_{\mathrm{o}}^{\*2}} - \sqrt{\left( \frac{\mathbf{R}}{\mathrm{H}} \right)^{2} + \mathrm{h}\_{\mathrm{s}}^{\*2}}}{\left( \mathrm{h}\_{\mathrm{c}}^{\*} - \mathrm{h}\_{\mathrm{o}}^{\*} \right)^{2} - \left( \mathrm{h}\_{\mathrm{c}}^{\*} - \mathrm{h}\_{\mathrm{s}}^{\*} \right)^{2}}. \tag{44}$$

Glover formula is derived from Equation (44) by making h∗ <sup>c</sup> = 1, h<sup>∗</sup> <sup>o</sup> = 0yh<sup>∗</sup> <sup>s</sup> = 1:

$$\mathbf{C} = \mathrm{asinh}\left(\frac{\mathrm{H}}{\mathrm{R}}\right) + \frac{\mathrm{R}}{\mathrm{H}} - \sqrt{\left(\frac{\mathrm{R}}{\mathrm{H}}\right)^{2} + 1}.\tag{45}$$

#### *3.3. The Reynolds and Elrick Model*

This model proposed by Reynolds and Elrick [27] assumes h∗ <sup>c</sup> = 1/2, h<sup>∗</sup> <sup>i</sup> = 0yh<sup>∗</sup> <sup>s</sup> = 1/2:

$$\mathbf{C} = 4 \left[ \frac{1}{2} \text{asinh} \left( \frac{\mathbf{H}}{2\mathbf{R}} \right) + \frac{\mathbf{R}}{\mathbf{H}} - \sqrt{\left( \frac{\mathbf{R}}{\mathbf{H}} \right)^2 + \frac{1}{4}} \right]. \tag{46}$$

#### *3.4. A Model for Stratified Porous Media*

Glover's model can be adapted for the case of stratified porous media. The well is considered to be in a medium composed of N layers of thickness Pj, j = 1,2, ... , N; the total hydraulic head, denoted as HT, is the height of the water column counted from the base of the well to the upper border of the N-th stratum.

The flow infiltrated by the walls of the j-th stratum is provided by Equation (43) modified as:

$$\mathbf{Q}\_{\rm oj} = \frac{2\pi \mathbf{K}\_{\rm s\parallel} \mathbf{p}\_{\rm j}^2}{\mathbf{C}\_{\rm j}},\tag{47}$$

where Ksj and Cj are the saturated hydraulic conductivity and the shape coefficient of the j-th stratum, respectively.

The shape coefficient is derived from Equation (44) denoting by Hj the hydraulic head at the base of the j-th stratum:

$$\mathbf{C}\_{\mathbf{j}} = \mathbf{h}\_{\mathbf{p}\mathbf{j}}^{\*} \frac{\mathbf{h}\_{\mathbf{C}\mathbf{j}}^{\*} \left[ \mathbf{as} \text{sinh} \left( \frac{\mathbf{p}\_{\mathbf{j}}}{\mathbf{K}} \mathbf{h}\_{\mathbf{s}\mathbf{j}}^{\*} \right) - \mathbf{as} \text{sinh} \left( \frac{\mathbf{p}\_{\mathbf{j}}}{\mathbf{K}} \mathbf{h}\_{\mathbf{c}\mathbf{j}}^{\*} \right) \right] + \sqrt{\left( \frac{\mathbf{R}}{\mathbf{P}\_{\mathbf{j}}} \right)^{2} + \mathbf{h}\_{\mathbf{c}\mathbf{j}}^{\*2}} - \sqrt{\left( \frac{\mathbf{R}}{\mathbf{P}\_{\mathbf{j}}} \right)^{2} + \mathbf{h}\_{\mathbf{s}\mathbf{j}}^{\*2}}}{\left( \mathbf{h}\_{\mathbf{c}\mathbf{j}}^{\*} - \mathbf{h}\_{\mathbf{c}\mathbf{j}}^{\*} \right)^{2} - \left( \mathbf{h}\_{\mathbf{c}\mathbf{j}}^{\*} - \mathbf{h}\_{\mathbf{s}\mathbf{j}}^{\*} \right)^{2}},\tag{48}$$

where h∗ Pj = Pj/Hj, h<sup>∗</sup> cj = hcj/Pj, h<sup>∗</sup> oj = hoj/Pj, h<sup>∗</sup> sj = hsj/Pj. It is noted that hcj, hoj y hsj are calculated from the base of the j-th stratum.

The Reynolds and Elrick model assumes h∗ <sup>c</sup> = 1/2, h<sup>∗</sup> <sup>o</sup> = 0yh<sup>∗</sup> <sup>s</sup> = 1/2 and therefore:

$$\mathbf{C}\_{\circ} = 4 \left[ \frac{1}{2} \text{asinh} \left( \frac{\mathbf{P}\_{\circ}}{2\mathbf{R}} \right) + \frac{\mathbf{R}}{\mathbf{P}\_{\circ}} - \sqrt{\left( \frac{\mathbf{R}}{\mathbf{P}\_{\circ}} \right)^{2} + \frac{1}{4}} \right] \mathbf{h}\_{\text{P}\text{\%}}.\tag{49}$$

The total flow is obtained as:

$$\mathbf{Q} = \sum\_{j=1}^{N} \mathbf{Q}\_{\mathbf{o}j} + \pi \mathbf{R}^2 \mathbf{K}\_{\mathbf{s1}} \tag{50}$$

where the flow at the bottom of the well has been added.

#### *3.5. Aplications*

To show the versatility of the solution, data obtained from an infiltration well built on the Queretaro Valley aquifer of radio are used: R = 0.3937 m (15.5) and depth PT = 36 m; Five strata were located in the profile (Figure 1). As drilling was carried out, the infiltration tests per stratum were carried out until the permanent regime was reached. The measured data are concentrated in Table 1, and the saturated hydraulic conductivity calculated from Equation (43) is also shown.

**Table 1.** Calculation of hydraulic conductivity per stratum, Equation (43).


Table 2 shows the flow rates calculated for each stratum when the well is full. In the last row is the total flow, Equation (50), and the saturated hydraulic conductivity corresponding to an equivalent homogeneous stratum.

**Figure 1.** Simplified well scheme.


**Table 2.** Calculation of the flow per stratum corresponding to a full well, Equations (47) and (49).

The artificial volume that can be recharged to the aquifer is 6972 L/s (602.3808 m3/day), however, stratum 1 only contributes 1.22% of the entire volume and is the deepest layer of the entire well (12 m). With these data, the cost benefit is analyzed and the decision is made to drill wells to a depth of 24 m with the understanding that we would contribute only 595.0318 m3/day to the aquifer but we reduce time and money.

The time to drill the well up to 36 m is 27 days at a cost of 25,457 USD. Therefore, if you choose to build two, the time taken would be 60 days with a total of 50,915 USD. Conversely, if you only drill to a depth of 24 m, the time taken is 7 days and a cost of 13,376 USD, which gives us a total cost of 53,504 USD for four wells drilled in 30 days. This is due to the fact that the material in the last 12 m is basalt and drilling progress is slower.

Regarding the volumes of recharge to the aquifer, with the four wells in the same area it would be 2380.1272 m3/day compared to 1204.7616 m3/day that we would obtain with only two at a depth of 36 m. Finally, the cost-benefit of annual recharge in the aquifer would be 16.24 m3/USD invested in four wells, compared to the 8.64 m3/USD that you have if you choose two wells.

#### **4. Conclusions**

In recent years the construction of artificial wells to recharge aquifers has been very popular in Mexico, however, as has been demonstrated in this work, the construction of a well at a greater depth does not necessarily give us a greater volume of recharge. The lack of information to calculate the total volumes has led to decisions being made with unscientific bases and, on several occasions, it has resulted in not achieving the expectations for which they were built.

This work provides a tool for knowing the behavior of the water infiltration rate in the porous medium in stratified media to artificially recharge an aquifer through wells. The analysis takes into account all the characteristics of the soil profiles that construct the well, resulting in the cost-benefit analysis of the complete operation to make a better decision.

It is widely demonstrated in the literature that several experimental tests are needed to know the behavior of this phenomenon, and that in order to know the process in detail, other factors that are not analyzed here must be taken into account: preferential flow in heterogeneous soils, trapped air, sediment deposit, among others.

When exploration drilling is done to propose a series of wells to recharge the aquifer, the data from the first layer is usually measured to simulate the behavior of the entire profile as a homogeneous stratum. However, as verified in this work, it is necessary to know the behavior of the entire well by stratum so that pertinent decisions are made, since a deeper well does not necessarily imply a greater volume of recharge to the aquifer.

**Author Contributions:** Conceptualization, C.F.; methodology, C.F.; software, C.F., A.Q., J.T.-A. and S.F.; validation, C.F., C.C. and A.Q.; formal analysis, C.F., C.C. and A.Q.; investigation, C.F.; data curation, J.T.-A. and S.F.; writing—original draft preparation, C.F.; writing—review and editing, C.C. and J.T.-A. All authors have read and agreed to the published version of the manuscript

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

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