Fuzzy solution of nonlinear Boussinesq equation

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.

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₀. Let and with ⁠, both differentiable at x₀. According to Bede & Stefanini (2013), it is:

• is (i)-gH-differentiable at x₀ if

  • (i)

    

• is (ii)-gH-differentiable at x₀ if

  • (ii)

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₀, t₀) 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₀, t₀) with respect to x. It is valid (Khastan & Nieto 2010; Allahviranloo et al. 2015) that:

• is [(i)-p]-differentiable with respect to x if:

  
• is [(ii)-p]-differentiable with respect to x if:

  

The second condition implies

Estimation of the term ⁠.

Table 1 shows the values of F_c, as well those of for different values of ⁠.

Table 1

Values of F_c and for different values of

Μ .	Fc .	.
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

Μ .	Fc .	.
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

View Large

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.
  The interval is divided into N equal parts:

  
• 2.

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

• 3.
  As ⁠, the value of arises from relation (24):

  
• 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_N = 1, where 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.
  A termination criterion is chosen. The following residual is formed:

  

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

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:

Remark

• For α = 0, the boundary condition is

• For α = 0, the boundary condition is

• For α = 1, the boundary condition is

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:

Table 4

Grid of N = 20 nodes for every case

.	.	.	.	.	.	.	.	.	.	.	.
	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

View Large

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.

Table 5

Different form values of Lockington (L), this study(ts) and Tolikas et al. (T)

.	.	μ .	.	.	.
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

View Large

Table 6

Relative errors

.	.
0.00477	0.000353
0.000737	0.000126
0.000836	0.000112

.	.
0.00477	0.000353
0.000737	0.000126
0.000836	0.000112

View Large

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,¹⁰–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.

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 − α.

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.

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

The authors declare there is no conflict.