## Abstract

In this paper, the solution of the one-dimensional second-order unsteady nonlinear fuzzy partial differential Boussinesq equation is examined. The physical problem concerns unsteady flow in a semi-infinite unconfined aquifer bordering a lake. In the examined problem, there is a sudden rise and subsequent stabilization of the lake's water level, thus the aquifer is recharging from the lake. The aquifer boundary conditions are considered fuzzy and, therefore, ambiguities are created to the solution of the overall physical problem. Then, the fuzzy problem is translated to a system of crisp boundary value problems. By using a Boltzmann transformation, the crisp problem is transformed into an integro-differential equation and solved with the help of a special numerical method. This method has a simple iterative procedure, which converges rapidly and is proven very accurate in comparison with other analytical methods. Additionally, the algebraic equation estimates very close to the storage coefficient, the flux at the stream-aquifer origin, and the water stored volume with respect to other exact solutions. This method is very useful for practical cases (artificial recharge problems, irrigation, and drainage projects), giving the possibility to decision-makers to take the right decision.

## HIGHLIGHTS

Solution of the fuzzy unsteady nonlinear flow.

Application of recharging an aquifer.

Numerical solution of the fuzzy integro-differential equation.

Accurate estimation of recharge coefficient.

Accurate estimation of stored volume/width and flux boundary at

*x*= 0.Application of a special numerical algorithm.

## INTRODUCTION

This equation was presented by Boussinesq (1904) with the following assumptions: (a) the inertial forces are negligible and (b) the horizontal component of velocity *u _{x}* does not vary depending on the depth and is a function of

*x*and

*t*. In 1904, Boussinesq provided a special solution to this nonlinear equation in the French journal ‘Journal de Mathématiques Pures et Appliquées’. Boussinesq's exact solution concerned the case of an aquifer overlying an impermeable layer, with boundary conditions similar to those of soil drained by a drain installed in the impermeable substratum. Polubarinova Kochina (1948, 1949, 2015) published a solution to Boussinesq's equation using the method of small disturbances. Tolikas

*et al.*(1984) obtained an approximate closed-form solution by applying similarity transformation and polynomial approximation. Lockington (1997) provided a simple approximate analytical solution using a weighted residual method. This method was applied to both the recharging and discharging of an unconfined aquifer due to a sudden change in the head at the origin. Lockington

*et al.*(2000) analyzed the problem of flow in a one-dimensional semi-infinite horizontal aquifer, initially dry and a head expressed as a power-law function of time at the origin. A quadratic approximate solution was derived in agreement with numerical results. Moutsopoulos (2010) applied Adomian's decomposition method and obtained a simple series solution with a few terms, and performing a benchmark test showed the advantages of his solution. Basha (2013) used the traveling wave method to obtain a nonlinear solution of a simple logarithmic form. The solution is adaptable to any flow situation that is recharge or discharge and allows the results of practical importance in

*Hydrology*. Additionally, algebraic equations are included for the velocity of the propagation front, wetting front position, and relationship for aquifer parameters. Chor

*et al.*(2013) provided a series solution for the nonlinear Boussinesq equation in terms of the Boltzmann transform in a semi-infinite domain. Jiang & Tang (2015) used a decomposition method, separating the original problem into linear diffusion equations. They developed a general approximate explicit solution in terms of error functions. More recently, Hayek (2019) provided an approximate solution, by introducing an empirical function with four parameters. The parameters were obtained using a numerical fitting procedure performed with Microsoft Excell Solver. Furthermore, analytical approaches are based mainly on Caputo fractional derivatives explored by authors to provide solutions to nonlinear partial differential equations. Specifically, Khan

*et al.*(2019) proposed a hybrid methodology of Shehu transformation along with the Adomian decomposition method while Shah

*et al.*(2019) and Rashid

*et al.*(2021) provide solutions to a system of nonlinear fractional Kortweg-de Vries partial differential equations based on Caputo operator, Shehu decomposition method, and the Shehu iterative transform method. Iqbal

*et al.*(2022) used a novel iterative transformation technique and homotopy perturbation transformation technique to calculate the fractional-order gas dynamics equation. Tzimopoulos

*et al.*(2021) used a transformed method of Wiedeburg (1890) to solve the one-dimensional Boussinesq equation for both the recharging and discharging of a homogeneous unconfined aquifer. Several other authors including Chen

*et al.*(1995), Parlange

*et al.*(2000), Telyakovskiy

*et al.*(2002), Pistiner (2008), Olsen & Telyakovskiy (2013), Bartlett & Porporato (2018), Tzimopoulos

*et al.*(2019) provide useful insight into the solution and are valuable tools for testing the accuracy of numerical methods.

Many numerical solutions to the problem of the water response to recharge or discharge of an aquifer have been developed in the past. Remson *et al.* (1971) give detailed information about the existing numerical methods for solving problems in subsurface hydrology. Tzimopoulos & Terzides (1975) investigated the case of water movement through soils drained by parallel ditches. Numerical solutions are presented based on implicit computational schemes of the Crank–Nicolson, Laasonen, and Douglas types. Experimental data, obtained by a Hele-Shaw model in the Laboratory, are in very good agreement with the values computed numerically by these implicit schemes. Tzimopoulos & Terzides (1976) studied a physical problem of a free surface flow toward a ditch or a river. The solution to Boussinesq's equation was made by the application of the finite-element method. Galerkin's method was used, leading to a system of nonlinear equations. Numerical results were compared with the exact solution of Boussinesq and with other finite-difference approximations. Frangakis & Tzimopoulos (1979) investigated a numerical model based on Boussinesq's equation describing the unsteady groundwater flow on impervious sloping bedrock. The numerical model uses the finite-element technique with Galerkin's method. The stability and accuracy of the method have been proved by the comparison of numerical results with the Crank–Nicolson scheme. Tzimopoulos & Tolikas (1980) investigated the problem of artificial groundwater recharge in the case of an unconfined aquifer described by the Boussinesq equation. The problem was solved by analytical and numerical methods. The finite-element method was used with square elements. Chávez *et al.* (2011) considered a problem of agricultural drainage described by the Boussinesq equation. They implemented an implicit numerical scheme in which interpolation parameters were used, that is, for space and time, respectively. Two discretization schemes of the time derivative were found: the mixed scheme and the head scheme. Both schemes were validated with one analytical solution. Bansal (2012, 2016) investigated the case of groundwater fluctuations in sloping aquifers induced by replenishment and seepage from a stream. For this case, the Boussinesq equation has been discretized using the MacCormack scheme. This scheme is a predictor–corrector scheme in which the predicted value of the head is obtained by replacing the spatial and temporal derivatives with forward differences. The corrector scheme is obtained by replacing the spatial derivative with a backward difference and the time derivative with a forward difference. Borana *et al.* (2013) employed a Crank–Nicolson finite-difference scheme to solve the Boussinesq equation for the case of infiltration phenomenon in a porous medium. They concluded that the Crank–Nicolson scheme is consistent with the physical phenomenon, and it is stable without having any restrictions on the stability ratio. Bansal (2017) investigated the case of the interaction of surface and groundwater in a stream-aquifer system. He derived a new analytical solution and employed the Du Fort and Frankel scheme for the comparison. The proposed scheme is an explicit finite-difference numerical scheme, proceeding in three-time levels. Nguyen (2018) investigated the case of stratified heterogeneous porous media. The model is a system of two equations: one for the water level in fissured porous blocks and one for the water level in system cracks. The discretized schemes are explicit in both cases. More recently, Samarinas *et al.* (2018, 2021, 2022) developed the Crank–Nicolson scheme in a fuzzy environment and solves the linearized fuzzy Boussinesq equation.

Since the aforementioned problem concerns differential equations, which present particular problems regarding fuzzy logic, a significant number of research studies were carried out in that field, especially regarding the fuzzy differentiation of functions. Initially, fuzzy differentiable functions were studied by Puri & Ralescu (1983), who generalized and extended Hukuhara's fundamental study (Hukuhara 1967) (H-derivative) of a set of values appearing in fuzzy sets. Kaleva (1987); Seikkala (1987) developed a theory on fuzzy differential equations. In the last years, several studies have been carried out in the theoretical and applied research field on fuzzy differential equations with an H-derivative (Kaleva 1987; Vorobiev & Seikkala 2002; O'Regan *et al.* 2003; Nieto & Rodríguez-López 2006). Nevertheless, in many cases, this method has presented certain drawbacks, since it has led to solutions with increasing support, along with increasing time (Diamond 2002). This proves that, in some cases, this solution is not a representative generalization of the classic case. To overcome this drawback, the generalized derivative gH (gH-derivative) was introduced (Bede & Gal 2005; Bede 2006; Stefanini 2010; Allahviranloo *et al.* 2015). The gH-derivative will be used henceforth, for a more extensive degree of fuzzy functions than the Hukuhara derivative.

In the present article, the problem of *fuzzy unsteady nonlinear* flow in a semi-infinite unconfined aquifer bordering a lake is introduced. In the examined problem, there is a sudden rise and subsequent stabilization of the lake's water level, thus the aquifer is recharging from the lake. The aquifer boundary conditions are considered fuzzy and therefore ambiguities are created in the solution of the problem. The fuzzy problem can be translated to a system of crisp boundary value problems hereafter called the corresponding system for the fuzzy problem. Subsequently, the crisp problem is solved with a numerical method based on Philip's method (1957), in which the initial Boussinesq equation is transformed with the help of a Boltzmann transformation into an integro-differential equation. The subsequent equation can readily be solved by changing the axis and by means of forward integration with one initial estimation of the stored volume determined by trial and improved by iteration. The extremely simple iterative procedure converges rapidly as it is explained in the Application section. This method proved to be very accurate in comparison with other analytical methods and also with the Runge–Kutta numerical method. Additionally, the algebraic equation estimates very close to the storage coefficient, the flux at the stream-aquifer boundary *x* = 0, and the water stored volume with respect to other exact solutions.

## METHODS

### Crisp model

### Dimensionless variables

## FUZZY MODEL DEFINITIONS

### Definitions

*Notation:* In order to facilitate the reader's nonfamiliarity with the fuzzy theory, some definitions are described below, concerning some preliminaries in fuzzy theory (Zadeh 1965) and some definitions about the differentiability.

*Definition 1.* A fuzzy set on a universe set *X* is a mapping , assigning to each element a degree of membership The membership function is also defined as with the properties:

(i) is upper semi-continuous, (ii) = 0, outside of some interval [c, d], (iii) there are real numbers , such that is increasing on [c, a], decreasing on [b, d] and = 1 for each (iv) is a convex fuzzy set (i.e., .

*Definition 2*. Let *X* be a Banach space and be a fuzzy set on *X*. The *α*-cuts of as , are defined and for *α* = 0 the closure is: .

*Definition 3*. Let (*X*) be the family of all nonempty compact convex subsets of a Banach space. A fuzzy set on *X* is called compact if (*X*), The space of all compact and convex fuzzy sets on *X* is denoted as (*X*).

*Definition 4.* Let . The *α*-cuts of , are: . According to the representation theorem of Negoita & Ralescu (1975) and the theorem of Goetschel & Voxman (1986), the membership function and the *α*-cut form of a fuzzy number are equivalent and in particular the *α*-cuts uniquely represent , provided that the two functions are monotonic (increasing, decreasing) and

*Definition 5*. gH- Differentiability (Bede & Stefanini 2013). Let be such that . Suppose that the functions , are real-valued functions, differentiable with respect to *x*, uniform with respect to. Then the function is gH-differentiable at a fixed if and only if one of the following two cases holds:

- (a)
is increasing, is decreasing as functions of

*α*, and , or - (b)
is decreasing, is increasing as functions of

*α*, and .

*Notation*: . In both of the above cases, the derivative is a fuzzy number.

*Definition 6*. gH-differentiable at *x*_{0}. Let and with , both differentiable at *x*_{0}. According to Bede & Stefanini (2013), it is:

is (i)-gH-differentiable at

*x*_{0}if- (i)

- (i)
is (ii)-gH-differentiable at

*x*_{0}if(ii)

*Definition 7*. g-differentiability. Let be such that . If and are differentiable real-valued functions with respect to

*x*, uniform for , then is g-differentiable and it is valid (Bede & Stefanini 2013) that:

*Definition 8.*The gH-differentiability implies g-differentiability, but the inverse is not true.

*Definition 9*. [gH-p] differentiability. A fuzzy-valued function of two variables is a rule that assigns to each ordered pair of real numbers (*x*, *t*) in a set , a unique fuzzy number denoted by . Let : (x_{0}, t_{0}) and , are real-valued functions and partial differentiable with respect to *x*. According to Khastan & Nieto (2010), Allahviranloo *et al.* (2015), it is:

*Definition 10*. Let : (), and be [gH-*p*]-differentiable at (*x*_{0}, *t*_{0}) with respect to *x*. It is valid (Khastan & Nieto 2010; Allahviranloo *et al.* 2015) that:

### Fuzzy nonlinear Boussinesq model

*et al.*(2015), and Tzimopoulos

*et al.*(2019) translating the above fuzzy problem to a system of the second order of crisp boundary value problems, hereafter called the

*corresponding system*for the fuzzy problem. Therefore, eight crisp boundary value problems (BVPs) systems are possible for the fuzzy problem {(1, 1), (1,2), (1,3), (1,4), (2, 1), (2,2), (2,3), (2,4)} with the same initial and boundary conditions.

*Note:*The problem is hereby restricted to the solution of the (1,1) system, which is described in detail. Case (1.1) is applied, since it gives a physical solution to the problem of aquifer recharging from the lake. Figure 1 presents the case of a lake recharged by an aquifer and Figure 2 presents the membership function of the boundary condition H(0,t) with r = 0.15.

### Solution of system (1,1)

#### First case

The second condition implies

Estimation of the term .

*F*is a weighted average:

_{c}Table 1 shows the values of *F _{c}*, as well those of for different values of .

Μ
. | F
. _{c} | . |
---|---|---|

1.25 | 1.17 | 1.16 |

1.5 | 1.33 | 1.32 |

2 | 1.67 | 1.63 |

3 | 2.33 | 2.26 |

4 | 3.00 | 2.89 |

5 | 3.67 | 3.52 |

Μ
. | F
. _{c} | . |
---|---|---|

1.25 | 1.17 | 1.16 |

1.5 | 1.33 | 1.32 |

2 | 1.67 | 1.63 |

3 | 2.33 | 2.26 |

4 | 3.00 | 2.89 |

5 | 3.67 | 3.52 |

### Stored volume

*et al.*(1971), the integral is equal to:and the stored volume can now be deduced from the following relation:

### Solving the integro-differential equation

For the convenience of the reader, the process for solving the above integro-differential equation is presented with the following flowchart in the explanatory form:

- 1.
- 2.
An initial estimation for value is assumed (Equations (25) and (30)).

- 3.
- 4.
Equation (21) is then used to compute and the process is continued.

- 5.
The computed value is compared now with the value . This last value has been computed from a nearly exact solution of Equation (10) in the vicinity of

*F*= 1, where_{N}*F*is considered a constant (for see Supplementary Material, appendix). - 6.
If does not agree closely enough with , a new value is chosen for and the procedure is repeated. The terme is (see Appendix): where

*A*(*w*) function is given in Table 2 and presented in Figure 4. - 7.

w
. | A(w)w
. | W
. | A(w)
. | w
. | A(w)
. |
---|---|---|---|---|---|

0 | 0 | 1.5 | 0.763 | 5 | 0.964 |

0.2 | 0.199 | 2 | 0.836 | 6 | 0.974 |

0.4 | 0.353 | 2.5 | 0.881 | 7 | 0.981 |

0.6 | 0.472 | 3 | 0.911 | 8 | 0.985 |

0.8 | 0.566 | 3.5 | 0.930 | 9 | 0.988 |

1 | 0.640 | 4 | 0.946 | 10 | 0.990 |

w
. | A(w)w
. | W
. | A(w)
. | w
. | A(w)
. |
---|---|---|---|---|---|

0 | 0 | 1.5 | 0.763 | 5 | 0.964 |

0.2 | 0.199 | 2 | 0.836 | 6 | 0.974 |

0.4 | 0.353 | 2.5 | 0.881 | 7 | 0.981 |

0.6 | 0.472 | 3 | 0.911 | 8 | 0.985 |

0.8 | 0.566 | 3.5 | 0.930 | 9 | 0.988 |

1 | 0.640 | 4 | 0.946 | 10 | 0.990 |

#### Second case

The second condition also implies The further solution is similar to Case 1.

## APPLICATION

In order to test the efficiency of the present study with other methods, the Runge–Kutta method (Tzimopoulos *et al.* 2021) and Lockington (1997) solutions are used. The last one has been fuzzified as follows:

### Crisp Lockington solution

### Fuzzy Lockington solution

Remark

For

*α*= 0, the boundary condition isFor

*α*= 0, the boundary condition isFor

*α*= 1, the boundary condition is

*K*= 20.00 m/d,

*S*= 0.27,

*h*

_{1}= 3 m,

*h*

_{0}= 2 m,

*t*= 5 d, and

*t*= 20d. In the special case where

*r*= 0.15 m and

*h*

_{0}= 2 m, , , , the values of ,, have been calculated from Lockington's solution and are presented in Table 3:

I . | 0 . | 1 . | 2 . |
---|---|---|---|

5.3453 | 5.4815 | 5.2837 | |

0.5215 | 0.4612 | 0.5749 | |

3.0000 | 2.5500 | 3.4500 |

I . | 0 . | 1 . | 2 . |
---|---|---|---|

5.3453 | 5.4815 | 5.2837 | |

0.5215 | 0.4612 | 0.5749 | |

3.0000 | 2.5500 | 3.4500 |

For every case, a grid of *N* = 20 nodes is used, which is presented in Table 4. All computations have been performed with a FORTRAN program, and the number of iterations was 7–8:

. | . | . | . | . | . | . | . | . | . | . | . |
---|---|---|---|---|---|---|---|---|---|---|---|

2,550 | 2,523 | 2,495 | 2,468 | 2,440 | 2,413 | 2,385 | 2,358 | 2,330 | 2,303 | 2,275 | |

3,000 | 2,950 | 2,900 | 2,850 | 2,800 | 2,750 | 2,700 | 2,650 | 2,600 | 2,550 | 2,500 | |

3,450 | 3,378 | 3,305 | 3,233 | 3,160 | 3,088 | 3,015 | 2,943 | 2,870 | 2,798 | 2,725 | |

. | . | . | . | . | . | . | . | . | . | . | . |

2,4550 | 2,400 | 2,350 | 2,300 | 2,250 | 2,200 | 2,150 | 2,100 | 2,050 | 2,000 | ||

2,248 | 2,220 | 2,193 | 2,165 | 2,138 | 2,110 | 2,083 | 2,055 | 2,028 | 2,000 | ||

2,653 | 2,580 | 2,508 | 2,435 | 2,363 | 2,290 | 2,218 | 2,145 | 2,073 | 2,000 |

. | . | . | . | . | . | . | . | . | . | . | . |
---|---|---|---|---|---|---|---|---|---|---|---|

2,550 | 2,523 | 2,495 | 2,468 | 2,440 | 2,413 | 2,385 | 2,358 | 2,330 | 2,303 | 2,275 | |

3,000 | 2,950 | 2,900 | 2,850 | 2,800 | 2,750 | 2,700 | 2,650 | 2,600 | 2,550 | 2,500 | |

3,450 | 3,378 | 3,305 | 3,233 | 3,160 | 3,088 | 3,015 | 2,943 | 2,870 | 2,798 | 2,725 | |

. | . | . | . | . | . | . | . | . | . | . | . |

2,4550 | 2,400 | 2,350 | 2,300 | 2,250 | 2,200 | 2,150 | 2,100 | 2,050 | 2,000 | ||

2,248 | 2,220 | 2,193 | 2,165 | 2,138 | 2,110 | 2,083 | 2,055 | 2,028 | 2,000 | ||

2,653 | 2,580 | 2,508 | 2,435 | 2,363 | 2,290 | 2,218 | 2,145 | 2,073 | 2,000 |

## RESULTS AND DISCUSSION

### Recharge coefficient-stored volume

*et al.*(1990), an accurate form of stored volume for the Boussinesq equation is presented as follows:

*C*is the dimensionless storage coefficient,

_{r}^{T}*K*= the hydraulic conductivity(m/d),

*h*

_{0}(m) = the initial water head, and

*S*= the dimensionless storability,

*κ*= a function of water level at origin, and . These values are taken here as reference values, and consequently, all the respective forms will be transformed into a similar form. Lockington (1997) form will be presented as follows:

in which are parameters, , *C _{r}* is a dimensional storage coefficient, and is a new dimensionless storage coefficient.

Table 5 presents the different form values of Lockington (*L*), this study (ts), and Tolikas *et al.* (*T*). Also, Table 6 presents the relative errors , (*i* = *L*, ts). The values of this study are very close to those of Tolikas *et al.* (1990). According to Table 6, all methods give small relative errors.

. | . | μ
. | . | . | . |
---|---|---|---|---|---|

2 | 1.275 | 0.3378 | 0.3363 | 0.3362 | |

2 | 1.500 | 0.6478 | 0.6474 | 0.6473 | |

2 | 1.725 | 0.9887 | 0.9880 | 0.9879 |

. | . | μ
. | . | . | . |
---|---|---|---|---|---|

2 | 1.275 | 0.3378 | 0.3363 | 0.3362 | |

2 | 1.500 | 0.6478 | 0.6474 | 0.6473 | |

2 | 1.725 | 0.9887 | 0.9880 | 0.9879 |

. | . |
---|---|

0.00477 | 0.000353 |

0.000737 | 0.000126 |

0.000836 | 0.000112 |

. | . |
---|---|

0.00477 | 0.000353 |

0.000737 | 0.000126 |

0.000836 | 0.000112 |

### Profiles-stored volume differences

Figures 5 and 6 illustrate the profiles of this study vs. the profiles of Lockington's and Runge–Kutta heads for *t* = 5d and 20d and Figures 7 and 8 illustrate the membership functions of for *t* = 5d and 20d. It is shown that there are small differences, between the profiles of this study and Lockington's profiles, and at the same time, the mean relative error of the corresponding profiles is greater than the mean relative error of the stored volume. This phenomenon can be explained by the fact that Lockington's solution gives smaller water head fronts than those of this study, and for that reason, the respective head profiles present some relative errors vs. this study in order to conserve the mass balance. The differences between this study and the Runge–Kutta method are very small and it is impossible to be visible graphically.

Finally, Figures 9,^{10}–11 illustrate the stored volume and the flux flow vs. Lockington's method as well as the membership functions of this study for *t* = 5d, 20d, and 40d. As the differences between the two methods regarding the stored volume are very small, it is impossible graphically to be visible.

### Decision makers

Figure 9 gives the possibility to decision makers, to take the right decision for irrigation or water resources problems. For example, for *t* = 40d the stored volume has the values: for *α* = 0, [51, 153], for *α* = 1(99.7). For different values of *α* according to the possibility theory (Zadeh 1965; Dubois *et al.* 2004; Mylonas 2022), it is easy for the decision makers to choose the appropriate value of *α* between 0 and 1 and to take the corresponding value of stored volume with probability *P* greater than 1 − *α*.

## CONCLUSIONS

The Hukuhara (1967) theory with the generalized Hukuhara (gH-derivative), as well as its extension by Allahviranloo *et al.* (2015) to partial differential equations, allows researchers to solve practical problems, useful in engineering. It is now possible for engineers to take the fuzziness of various sizes involved into consideration when calculating and making the right decision on their work.

Boussinesq's nonlinear equation, regarding the aquifer case study, has a new fuzzy numerical solution, which is very accurate in comparison to the Runge–Kutta method and extremely easy and simple to calculate. The estimated storage coefficient and stored volume are very accurate in comparison to other methods.

The fuzzy water volume variations (spreads), regarding their average value, amount to the same quantity of the average values of all times and, therefore, equal to the initial relative fuzziness in boundary values. The fuzzy water flux flow tends asymptotically to zero, while the recharging fuzzy water volume tends asymptotically to infinity.

It is important here to note that for practical cases (artificial recharge problems, irrigation, and drainage projects), the engineers should take the right decision, knowing the deviations of the crisp value of water volume from the fuzzy ones, with a certain probability.

## DATA AVAILABILITY STATEMENT

All relevant data are included in the paper or its Supplementary Information.

## CONFLICT OF INTEREST

The authors declare there is no conflict.

## REFERENCES

*Applications in Fuzzy Statistic and Approximate Reasoning*

*Numerical and Analytical Analysis of Flow in Stratified Heterogeneous Porous Media*