**1. Introduction**

The horizontal water flow concerning unconfined aquifers without precipitation is described by the one-dimensional second-order unsteady nonlinear partial differential equation, called the Boussinesq equation:

$$\frac{\partial h}{\partial t} = \frac{K}{S} \frac{\partial}{\partial x} \left( h \frac{\partial h}{\partial x} \right) \tag{1}$$

where *K* = hydraulic conductivity (LT−1), *S* = effective porosity(L0T0), *h* = piezometric head (L), *x* = horizontal coordinate (L), *t* = time (T).

Boussinesq (1904) [1] first proposed the above equation according to the assumption that the horizontal component of velocity ux does not change with depth and is a function of x and t, while the inertial forces are negligible. A unique solution to this nonlinear equation was published by Boussinesq in the "*Journal de Mathématiques Pures et Appliquées*" in 1904. With boundary conditions such as those of soil drained by drains placed in the impermeable substratum, Boussinesq solution dealt with the case of an aquifer atop an impermeable layer. Using the small disturbance method, Polubarinova-Kochina (952, 1962) [2,3] published a

**Citation:** Tzimopoulos, C.; Samarinas, N.; Papadopoulos, K.; Evangelides, C. Fuzzy Analytical Solution for the Case of a Semi-Infinite Unconfined Aquifer. *Environ. Sci. Proc.* **2023**, *25*, 70. https://doi.org/ 10.3390/ECWS-7-14303

Academic Editor: Athanasios Loukas

Published: 3 April 2023

**Copyright:** © 2023 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 (https:// creativecommons.org/licenses/by/ 4.0/).

solution to Boussinesq's equation. By utilizing polynomial approximation and similarity transformation, Tolikas et al. (1984) [4] found an approximate solution. A weighted residual approach was used by Lockington (1997) [5] to provide an easily applicable analytical solution. Due to an abrupt change in the head at the origin, this approach was used for both the recharging and discharging of an unconfined aquifer. By using Adomian's decomposition method, Moutsopoulos (2010) [6] arrived at a simple series solution with a limited number of terms while he also performed a benchmark test, which demonstrated the benefits of his solution. The problem of flow in a one-dimensional semiinfinite horizontal aquifer with an initially dry head and a power-law function of time at the origin was examined by Lockington et al. in 2000 [7]. In accordance with the numerical outcomes, an approximative quadratic solution was developed. Using the traveling wave method, Basha (2013) [8] obtained a nonlinear solution with an easily applicable logarithmic form. The solution allows for the results to be of practical value in hydrology and is adaptable to any flow situation, whether it be recharge or discharge condition. There are also algebraic formulae for the propagation front velocity, the location of the wetting front, and the linkage between the characteristics of the aquifer. The nonlinear Boussinesq equation was given by a series solution by Chor et al. (2013) [9] in terms of the Boltzmann transform in a semi-infinite domain. An approximate solution was recently obtained by Hayek (2019) [10], who introduced an empirical function with four parameters. Using Microsoft Excel Solver, a numerical fitting approach was used to acquire the parameters. An approximate analytical solution for the recharge and discharge of a homogeneous unconfined aquifer was published by Tzimopoulos et al. in 2022 [11]. Numerous other studies [12–16] offer helpful clarification on the solution, offering a way for testing and accuracy of the numerical methods.

The definition of the initial flow condition, the method of linearizing the Boussinesq equation, the definition of drain spacing and hydraulic conductivity, boundary conditions, etc., are just a few examples of the ambiguities and uncertainties that the physical problem described by the Boussinesq equation presents [17,18]. Without taking into account the ambiguities and uncertainties of the groundwater flow problems, wrong managemen<sup>t</sup> decisions could be made, leading to a number of significant negative environmental, social, and economic effects. Fuzzy algorithms were used to solve this problem for all the aforementioned reasons.

The fuzzy logic theory is a useful tool for modeling ambiguity, developed by Lofti Zadeh (1965) [19]. Its development has had a significant impact on both theoretical problems [20–23] as well as engineering and hydraulic problems [24]. To solve fuzzy differential equations, some analytical and numerical approaches have recently been put forth. Chang and Zadeh (1972) [25] first introduced the concept of fuzzy derivative, while Dubois and Prade (1982) [26] followed by using the extension principle in their approach. Fuzzy differential functions were studied by Puri and Ralescu (1983) [27], who extended Hukuhara's derivative (H-derivative) [28] of a set of values appearing in fuzzy sets. Kaleva (1987, 1990) [29,30] and Seikkala (1987) [31] developed the fuzzy initial value problem, while Abbasbandy and Allahviranloo (2002) [32] presented a numerical algorithm for solving fuzzy ordinary differential equations based on the second Taylor method. However, this method presented certain drawbacks, and in many cases, this solution was not a good generalization of the classic case. The generalized Hukuhara differentiability (gH—differentiability) was introduced by [33,34], overcoming this drawback. This new derivative is defined for a larger class of fuzzy functions than the Hukuhara derivative. Allahviranloo et al. (2015) [35] introduced the (gH-p) differentiability for partial derivatives as an extension of the above theory. The gH-p differentiability was used by Tzimopoulos et al. (2018, 2020) [17,36], providing a fuzzy linear analytical solution to a parabolic partial differential equation, and also Tzimopoulos et al. (2018) [18] obtained a fuzzy linear analytical solution to the Boussinesq equation in the case of an unconfined aquifer problem.

In this paper, a comparison between the fuzzy analytical nonlinear Boussinesq equation and the Runge–Kutta numerical method is presented for the proposed fuzzy analytical

solution accuracy evaluation and effectiveness check. This comparison attests to the accuracy of the former. Additionally, an application of the Possibility Theory [37,38] to the *α*-cuts of the water level profiles, as well as of the water volume variation, provides the water levels and the water volume confidence intervals with probability *p=1* − *α*. Thus, a combination of fuzzy theory with the Possibility Theory allows managers and engineers to solve practical hydraulic problems, making the right decision.

#### **2. Constructing the Fuzzy Model and Solution**

### *2.1. Crisp Model*

As mentioned in the introduction to the case of one-dimensional horizontal flow, the Equation (1) corresponds to Boussinesq equation

The initial and boundaries conditions are as follows:

$$\begin{aligned} t &= 0, h(\mathbf{x}, 0) = h\_{0\prime} \\ t &> 0, h(0, t) = h\_{1\prime} \underset{\mathbf{x} \to \infty}{\text{h}}(\mathbf{x}, t) = h\_0 \end{aligned}$$

The solution [11] of the above Equation (1) is as follows:

$$h = h\_0 + (h\_0 - h\_1)\Omega(\mu, \xi).$$

where

$$\begin{array}{c} \Omega(\mu,\xi) = \ell F(\frac{x}{\xi}) - \Phi(\frac{x}{\xi}), \Phi(\frac{x}{\xi}) = \operatorname{erfc}(\frac{x}{\xi}),\\ \xi = \frac{x}{2\sqrt{\frac{h\_1}{\xi}}}, \ell = (h\_0 - h\_1)/h\_1, \mu = h\_1/h\_0, \\ F(\xi) = -\frac{1}{\pi} + \left(\frac{1}{2} + \frac{1}{\pi}\right)\Phi - \frac{1}{\sqrt{\pi}}(1 - \Phi)\xi e^{-\frac{\xi^2}{5}} + \frac{1}{\pi}(1 - e^{-2\xi^2}) - \frac{1}{2}\Phi^2. \end{array}$$

### *2.2. Fuzzy Model*

We write Equation (1), in its fuzzy form as follows:

$$\frac{\partial \tilde{h}}{\partial t} = \frac{\mathcal{K}}{\mathcal{S}} \frac{\partial}{\partial \mathbf{x}} (\tilde{h} \frac{\partial \tilde{H}}{\partial \mathbf{x}}) = \frac{\mathcal{K}}{\mathcal{S}} \{ (\frac{\partial \tilde{h}}{\partial \mathbf{x}})^2 + \tilde{h} \frac{\partial^2 \tilde{h}}{\partial \mathbf{x}^2} \}, \tag{2}$$

with the new boundary and initial conditions:

$$\begin{array}{ll}\text{Initial conditions} & \text{Boundary conditions} \\ t = 0, \check{h}(\mathbf{x}, 0) = \check{h}\_{0\prime} t > 0, \quad \check{h}(0, t) = \check{h}\_{1\prime} \underset{\mathbf{x} \to \infty}{\check{h}}(\mathbf{x}, t) = \check{h}\_{0} \end{array} \tag{3}$$

By converting the aforementioned fuzzy problem into a system of second order crisp boundary value problems—(referred as the *corresponding system* for the fuzzy problem)— and applying the theory of [18,35,39,40], we are able to solve the fuzzy problem (2), subject to the boundary and initial conditions (3). As a result, for the fuzzy problem, eight crisp BVPs systems are feasible under identical initial and boundary conditions.

Note: We shall now turn our attention to the first system's solution as it offers a practical solution to the case of the lake refilling the aquifer.

#### *2.3. Solution of the First System*

The nonlinear one-dimensional horizontal flow equations related to the first case (left boundary) and to the second case (right boundary) are provided with the following expressions:

$$\begin{array}{cccc} & \text{First case} & \text{Second case} \\ \frac{\partial h^{-}}{\partial t} = \frac{K}{\mathsf{S}} \left\{ h^{-} \frac{\partial^{2}h^{-}}{\partial x^{2}} + \left( \frac{\partial h^{-}}{\partial x} \right)^{2} \right\} & \frac{\partial h^{+}}{\partial t} = \frac{K}{\mathsf{S}} \left\{ h^{+} \frac{\partial^{2}h^{+}}{\partial x^{2}} + \left( \frac{\partial h^{+}}{\partial x} \right)^{2} \right\} \end{array} \tag{4}$$

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

We have used the data of Lockington (1997) [5], which means *K* = hydraulic conductivity = 20 m/d, *S* = effective porosity = 0.27, *h*1 = 3m, *h*0 = 2m. For the first case, the solution is as follows [11]:

$$h^-(\mu\_\prime \check{\mathfrak{z}}) = h\_0^- + (h\_0^- - h\_1^-) \Omega(\mu\_\prime \check{\mathfrak{z}})^2$$

Figure 1 illustrates the water level profiles for *t* = 5d of the new analytical method vs Runge–Kutta method (reference solution) and for *α* = 0, 1 value. The two methods approach each other closely. In addition, the water table in a period of 5 days approaches 125 m in length. If we consider also the hydraulic conductivity value (*K* = 20 m/d), the results are absolutely reasonable regarding the physical problem. Figure 2 illustrates the stored volume variation vs. time, and Figure 3 illustrates the membership function of *V* (*<sup>x</sup>*, *t*) for *t* = 5d. According to possibility theory [37,38], every function [*V* ] estimates the crisp function *V*, and the *α*-cut {*V* }*a* = [*<sup>V</sup>*<sup>−</sup>*a* , *<sup>V</sup>*+*a* ] should be interpreted as the confidence intervals of *V* with a probability *p* ≥ 1 − *a*. In this regard, in Figure 3, it is seen that for *α* = 0.05, the value of water volume lies in the interval [5.226, 14.310] with a probability higher than 95%, according to the possibility theory.

**Figure 1.** Water level profiles for *t* = 5d.

**Figure 2.** Stored volume variation vs. time.

**Figure 3.** Membership function of *V* (*<sup>x</sup>*, *t*).
