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.

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

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

This equation was presented by Boussinesq (1904) with the following assumptions: (a) the inertial forces are negligible and (b) the horizontal component of velocity ux 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.

Crisp model

As mentioned in the Introduction section, in the case of one-dimensional horizontal flow, the Boussinesq equation is
formula
(1)
The initial and boundary conditions are
formula
(2)

Dimensionless variables

We introduce now dimensional variables as follows:
formula
and the Boussinesq equation becomes:
formula
(3)
with new initial and boundary conditions:
formula
(4)

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

  • is (i)-gH-differentiable at x0 if

    • (i)

  • is (ii)-gH-differentiable at x0 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:
formula
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 : (x0, t0) and , are real-valued functions and partial differentiable with respect to x. According to Khastan & Nieto (2010), Allahviranloo et al. (2015), it is:

is [(i)-p]-differentiable with respect to x at (x0, t0) if:
formula
is [(ii)-p]-differentiable with respect to x at (x0, t0) if:
formula

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

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

Fuzzy nonlinear Boussinesq model

Εquation (3) is written in its fuzzy form as follows:
formula
(5)
with the new boundary and initial conditions:
formula
(6)
Solutions to the fuzzy problem (5) and the boundary and initial conditions (6) can be found utilizing the theory of Khastan & Nieto (2010), Bede & Stefanini (2013), Allahviranloo 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.
formula
formula
formula
formula
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.
Figure 1

A definition sketch of the investigated problem.

Figure 1

A definition sketch of the investigated problem.

Close modal
Figure 2

A membership function of .

Figure 2

A membership function of .

Close modal

Solution of system (1,1)

formula

First case

formula
(7)
Boundary conditions are as follows:
formula
(8)
Initial condition is as follows:
formula
(9)
In order to facilitate the calculations, it is posed that and the above equation becomes:
formula
(10)
Boundary conditions are as follows:
formula
(11)
Initial condition is as follows:
formula
(12)
Applying now the well-known Boltzmann similarity transformation:
formula
an ordinary differential equation (ODE) is obtained with new boundary conditions:
formula
(13)
formula
(14)

The second condition implies

According to this condition:
formula
and Equation (13) is transformed now to the following integro-differential equation:
formula
(15)
Equation (15) is referred now to as a new finite-difference system, presented in Figure 3, where the axis x is for function F and axis y is for :
Figure 3

A new finite-difference node system.

Figure 3

A new finite-difference node system.

Close modal
TermsTaking now a typical node r + 1/2, finite-difference approximations may be constructed for the form:
formula
(16)
A new term is now introduced:
formula
(17)
formula
(18)
The new term is denoted as follows:
formula
(19)
A similar form is constructed now for node r–1/2, and a second new form appears as follows:
formula
formula
(20)
The desired recurrence relation is now as follows:
formula
(21)
The second term of Equation (15) may be written according to Equation (18) as:
formula
(22)
which combined with Equation (19) and multiplied with ΔF i.e.,
formula
(23)
gives the new term of :
formula
(24)

Estimation of the term .

An estimation of this term is given by Kirkham & Powers (1972) as follows:
formula
(25)
where Fc is a weighted average:
formula
In this present case, this is equal to:
formula
(26)
Now the integral can be written as:
formula
(27)
because it is
formula
(28)
and the integral has its final form
formula
(29)

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

Table 1

Values of Fc and for different values of

ΜFc
1.25 1.17 1.16 
1.5 1.33 1.32 
1.67 1.63 
2.33 2.26 
3.00 2.89 
3.67 3.52 
ΜFc
1.25 1.17 1.16 
1.5 1.33 1.32 
1.67 1.63 
2.33 2.26 
3.00 2.89 
3.67 3.52 

The new values are calibrated in order to be in good agreement with the corresponding numerical values and the function has the following form:
formula
(30)

Stored volume

According to Remson et al. (1971), the integral is equal to:
formula
(31)
and the stored volume can now be deduced from the following relation:
formula
(32)
Finally, the stored volume is estimated as:
formula
(33)
and the flux boundary at x = 0 is equal to:
formula
(34)

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

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

  • 3.
    As , the value of arises from relation (24):
    formula
  • 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 FN = 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:
    formula

Table 2

Values of A(w) (Remson et al. 1971) and for w > 3

wA(w)wWA(w)wA(w)
1.5 0.763 0.964 
0.2 0.199 0.836 0.974 
0.4 0.353 2.5 0.881 0.981 
0.6 0.472 0.911 0.985 
0.8 0.566 3.5 0.930 0.988 
0.640 0.946 10 0.990 
wA(w)wWA(w)wA(w)
1.5 0.763 0.964 
0.2 0.199 0.836 0.974 
0.4 0.353 2.5 0.881 0.981 
0.6 0.472 0.911 0.985 
0.8 0.566 3.5 0.930 0.988 
0.640 0.946 10 0.990 
Figure 4

Function ?(w).

If the residual is sufficiently small i.e.,
formula
the procedure is terminated. Otherwise, the iterations are repeated, taking a new value for
formula
  • 8.

    In the end of the procedure, all the values of have been computed, corresponding to the values of

Second case

formula
(35)
Boundary conditions are as follows:
formula
(36)
Initial condition is as follows:
formula
(37)
In order to facilitate the computations, it is posed that and the abovementioned equation becomes:
formula
(38)
Boundary conditions are follows:
formula
(39)
The initial condition is as follows:
formula
(40)
Applying now the well-known Boltzmann similarity transformation:
formula
the following equation is obtained:
formula
(41)
formula
(42)

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:

Crisp Lockington solution

formula

Fuzzy Lockington solution

formula

Remark

  • For α = 0, the boundary condition is

  • For α = 0, the boundary condition is

  • For α = 1, the boundary condition is

Lockington's solution can now be written as follows:
formula
where has correspondingly the above values. The same example with Lockington (1997) and Moutsopoulos (2010) is used, namely the estimation of a water table for an aquifer with K = 20.00 m/d, S = 0.27, h1 = 3 m, h0 = 2 m, t = 5 d, and t = 20d. In the special case where r = 0.15 m and h0 = 2 m, , , , the values of ,, have been calculated from Lockington's solution and are presented in Table 3:
Table 3

Values of , and calculated from Lockington's solution

I012
 5.3453 5.4815 5.2837 
 0.5215 0.4612 0.5749 
 3.0000 2.5500 3.4500 
I012
 5.3453 5.4815 5.2837 
 0.5215 0.4612 0.5749 
 3.0000 2.5500 3.4500 
Flux Q at x = 0 and stored volume V of Lockington solution:
formula
formula

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  

The following Figures 5 and 6 illustrate the profiles of this study vs. Runge–Kutta and Lockington solutions for t = 5d and t = 20d.
Figure 5

Profiles of this study vs. Runge?Kutta and Lockington solutions for t = 5d.

Figure 5

Profiles of this study vs. Runge?Kutta and Lockington solutions for t = 5d.

Close modal
Figure 6

Profiles of this study vs. Runge–Kutta and Lockington solutions for t = 20d.

Figure 6

Profiles of this study vs. Runge–Kutta and Lockington solutions for t = 20d.

Close modal
The following Figures 7 and 8 illustrate the membership functions of this study vs. the Lockington solution for t = 5d and t = 20d.
Figure 7

Membership functions for t = 5d of this study vs. Lockington solution.

Figure 7

Membership functions for t = 5d of this study vs. Lockington solution.

Close modal
Figure 8

Membership functions of for t = 20d of this study vs. Lockington solution.

Figure 8

Membership functions of for t = 20d of this study vs. Lockington solution.

Close modal
Figures 911 illustrate the stored volume of this study vs. Lockington solution for different times, the membership functions of vs. Lockington and flux flow Q of this study vs. Lockington's solution, respectively.
Figure 9

Stored volume of this study vs. Lockington solution.

Figure 9

Stored volume of this study vs. Lockington solution.

Close modal
Figure 10

Membership functions of for t = 5d, 20d, and 40d.

Figure 10

Membership functions of for t = 5d, 20d, and 40d.

Close modal
Figure 11

The flux flow Q of this study vs. Lockington solution.

Figure 11

The flux flow Q of this study vs. Lockington solution.

Close modal

Recharge coefficient-stored volume

According to Tolikas et al. (1990), an accurate form of stored volume for the Boussinesq equation is presented as follows:
formula
in which CrT is the dimensionless storage coefficient, K = the hydraulic conductivity(m/d), h0(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:
formula
formula

in which are parameters, , Cr is a dimensional storage coefficient, and is a new dimensionless storage coefficient.

In our study, the stored volume is presented as follows (Equation (34)):
formula
formula
(30)

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)

μ
 1.275 0.3378 0.3363 0.3362 
 1.500 0.6478 0.6474 0.6473 
 1.725 0.9887 0.9880 0.9879 
μ
 1.275 0.3378 0.3363 0.3362 
 1.500 0.6478 0.6474 0.6473 
 1.725 0.9887 0.9880 0.9879 
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 

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

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.

Allahviranloo
T.
,
Gouyandeh
Z.
,
Armand
A.
&
Hasanoglu
A.
2015
On fuzzy solutions for heat equation based on generalized Hukuhara differentiability
.
Fuzzy Sets and Systems
265
,
1
23
.
https://doi.org/10.1016/j.fss.2014.11.009
Bansal
R. K.
2016
Approximate analytical solution of Boussinesq equation in homogeneous medium with leaky base
.
Applications and Applied Mathematics: An International Journal (AAM)
11
(
1
),
184
191
.
Bartlett
M. S.
&
Porporato
A.
2018
Class of exact solutions of the Boussinesq equation for horizontal and sloping aquifers
.
Water Resources Research
54
(
2
),
767
778
.
Basha
H. A.
2013
Traveling wave solution of the Boussinesq equation for groundwater flow in horizontal aquifers
.
Water Resources Research
49
,
1668
1679
.
doi:10.1002/wrcr.20168
.
Bede
B.
2006
A note on two-point boundary value problems associated with non-linear fuzzy differential equations
.
Fuzzy Sets and Systems
157
,
986
989
.
https://doi.org/10.1016/j.fss.2005.09.006
.
Bede
B.
&
Gal
S. G.
2005
Generalizations of the differentiability of fuzzy-number-valued functions with applications to fuzzy differential equations
.
Fuzzy Sets and Systems
151
,
581
599
.
https://doi.org/10.1016/j.fss.2004.08.001
Bede
B.
&
Stefanini
L.
2013
Generalized differentiability of fuzzy-valued functions
.
Fuzzy Sets and Systems
230
,
119
141
.
Borana
R. N.
,
Pradhan
V. H.
&
Mehta
N.
2013
Numerical solution of Boussinesq equation arising in one-dimensional infiltration phenomenon by using finite difference method
.
International Journal of Research in Engineering & Technology
2
(
8
),
202
209
.
Available from: http://www.ijret.org.
Boussinesq
J.
1904
Recherches théoriques sur l’écoulement des nappes d'eau infiltrées dans le sol et sur le débit des sources
.
Journal de Mathématiques Pures et Appliquées
10
(
5
),
5
78
.
Chávez
C.
,
Fuentes
C.
,
Zavala
M.
&
Brambila
F.
2011
Numerical solution of the Boussinesq equation. Application to the agricultural drainage
.
African Journal of Agricultural Research
6
(
18
),
4210
4222
.
Chen
Z. X.
,
Bodvarsson
G. S.
,
Witherspoon
E. A.
&
Yortsos
Y. C.
1995
An integral equation formulation for the unconfined flow of groundwater with variable inlet conditions
.
Transport in Porous Media
18
,
15
36
.
Chor
T.
,
Dias
N. L.
&
de Zarate
A. R.
2013
An exact series and improved numerical and approximate solutions for the Boussinesq equation
.
Water Resources Research
49
,
7389
7387
.
doi:10.1002/wrcr.20543,2013
.
Diamond
P.
2002
Brief note on the variation of constants formula for fuzzy differential equations
.
Fuzzy Sets and Systems
129
,
65
71
.
https://doi.org/10.1016/S01
Dubois
D.
,
Foulloy
L.
,
Mauris
G.
&
Prade
H.
2004
Probability possibility transformations, triangular fuzzy sets and probabilistic inequalities
.
Reliable Computing
10
,
273
297
.
Frangakis
C. N.
&
Tzimopoulos
C.
1979
Unsteady groundwater flow on sloping bedrock
.
Water Resources Research
15
,
176
180
.
https://doi.org/10.1029/WR015i001p00176
Goetschel
R.
&
Voxman
W.
1986
Elementary fuzzy calculus
.
Fuzzy Sets and Systems
18
,
31
43
.
Hukuhara
M.
1967
Intégration des applications mesurables dont la valeur est un compact convexe
.
Funkcialaj Ekvacioj
10
,
205
233
.
Iqbal
N.
,
Akgul
A.
,
Shah
R.
,
Bariq
A.
,
Al-Sawalha
M. M.
&
Ali
A.
2022
On solutions of fractional-order gas dynamics equation by effective techniques
.
Journal of Function Spaces
2022
,
1
14
.
https://doi.org/10.1155/2022/3341754
Kaleva
O.
1987
Fuzzy differential equations
.
Fuzzy Sets and Systems
24
,
301
317
.
doi:10.1016/0165-0114(87)90029-7
.
Khan
H.
,
Farooq
U.
,
Shah
R.
,
Baleanu
D.
,
Kumam
P.
&
Arif
M.
2019
Analytical solutions of (2 + Time fractional order) dimensional physical models, using modified decomposition method
.
Applied Sciences
10
(
1
),
1
20
.
Khastan
A.
&
Nieto
J. J.
2010
A boundary value problem for second order fuzzy differential equations
.
Nonlinear Analysis
72
,
3583
3593
.
Kirkham
D.
&
Powers
W. L.
1972
Advanced Soil Physics
.
Wiley-Interscience
,
New York
,
534
.
Lockington
D. A.
1997
Response of unconfined aquifer to sudden change in boundary head
.
Journal of Irrigation and Drainage Engineering
123
,
24
27
.
Lockington
D. A.
,
Parlange
J. Y.
,
Parlange
M. B.
&
Selker
J.
2000
Similarity solution of the Boussinesq equation
.
Advances in Water Resources
23
,
725
729
.
Mylonas
N.
2022
Applications in Fuzzy Statistic and Approximate Reasoning
.
PhD
,
Dimokritos University of Thrace-Greece
, p.
164
,
(In Greek)
.
Negoita
C. V.
&
Ralescu
D. A.
1975
Representation theorems for fuzzy concepts
.
Kybernetes
4
(
3
),
169
174
.
Nguyen
T.
2018
Numerical and Analytical Analysis of Flow in Stratified Heterogeneous Porous Media
.
M.Sc.
,
University of Stavanger
,
Norway
, p.
76
.
Nieto
J. J.
&
Rodríguez-López
R.
2006
Bounded solutions for fuzzy differential and integral equations
.
Chaos, Solitons and Fractals
27
,
1376
1386
.
doi:10.1016/j.chaos.2005.05.012
.
O'Regan
D.
,
Lakshmikantham
V.
&
Nieto
J. J.
2003
Initial and boundary value problem for fuzzy differential equations
.
Nonlinear Analysis
54
,
405
415
.
doi:10.1016/S0362-546X(03)00097-X
.
Olsen
J. S.
&
Telyakovskiy
A. S.
2013
Polynomial approximate solutions of a generalized Boussinesq equation
.
Water Resources Research
49
(
5
),
3049
3053
.
doi:10.1002/wrcr.20242,2013
.
Parlange
J. Y.
,
Hogarth
W. L.
,
Govindaraju
R. S.
,
Parlange
M. B.
&
Lockington
D.
2000
On an exact analytical solution of the Boussinesq equation
.
Transport in Porous Media
39
,
339
345
.
Philip
J. R.
1957
Numerical solution of equations of the diffusion type with diffusivity concentration-dependent. II
.
Australian Journal of Physics
30
,
20
42
.
Pistiner
A.
2008
Similarity solution to unconfined flow in an aquifer
.
Transport in Porous Media
71
,
265
272
.
doi:10.1007/s11242-007-9124-5
.
Polubarinova-Kochina
P. Y.
1948
On a non-linear differential equation occurring in seepage theory
.
Doklady Akadémii Nauk SSSR
35
(
6
).
Polubarinova-Kochina
P. Y.
1949
On unsteady motions of groundwater during seepage from water reservoirs
.
PMM (Prinkladaya Matematica I Mekhanica)
13
(
2
).
Polubarinova-Kochina
P. Y.
2015
Theory of Groundwater Movement, 2015 ed
, 1st edn.
Princeton University Press, NJ
,
1962, USA, 2015; pp.613; ISBN
.
Puri
M. L.
&
Ralescu
D. A.
1983
Differentials of fuzzy functions
.
Journal of Mathematical Analysis and Applications
91
,
552
558
.
Rashid
S.
,
Khalid
A.
,
Sultana
S.
,
Hammouch
Z.
,
Shah
R.
&
Alsharif
M. A.
2021
A novel analytical view of time-fractional korteweg-de Vries equations via a new integral transform
.
Symmetry
13
(
7
),
1
31
.
https://doi.org/10.3390/sym13071254
Remson
I.
,
Hornberger
G. M.
&
Moltz
F. G.
1971
Numerical Methods in Subsurface Hydrology
.
Wiley-Interscience
,
NY, London, Sydney, Toronto
, pp.
389
.
Samarinas
N.
,
Tzimopoulos
C.
&
Evangelides
C.
2018
Fuzzy numerical solution to horizontal infiltration
.
International Journal of Circuits, Systems and Signal Processing
12
,
326
332
.
Samarinas
N.
,
Tzimopoulos
C.
&
Evangelides
C.
2021
Fuzzy numerical solution to the unconfined aquifer problem under the Boussinesq equation
.
Water Supply
21
(
6
),
3210
3224
.
https://doi.org/10.2166/ws.2021.115
Samarinas
N.
,
Tzimopoulos
C.
&
Evangelides
C.
2022
An efficient method to solve the fuzzy Crank-Nicolson scheme with application to the groundwater flow problem
.
Journal of Hydroinformatics
24
(
3
),
590
609
.
https://doi.org/10.2166/hydro.2022.150
Seikkala
S.
1987
On the fuzzy initial problem
.
Fuzzy Sets and Systems
24
,
319
330
.
doi:10.1016/0165-0114(87)90030-3
.
Shah
R.
,
Khan
H.
,
Kuman
P.
&
Arif
M.
2019
An analytical technique to solve the system of nonlinear fractional partial differential equations
.
Mathematics
505
(
7
),
1
16
.
doi:10.3390/math7060505
.
Stefanini
L.
2010
A generalization of Hukuhara difference and division for interval and fuzzy Arithmetic
.
Fuzzy Sets and Systems
0
(
161
),
1564
1584
.
doi:10.1016/j.fss.2009.06.009
.
Telyakovskiy
A. S.
,
Braga
G. A.
&
Furtado
F.
2002
Approximate similarity solutions to the Boussinesq equation
.
Advances in Water Resources
25
(
2
),
191
194
.
doi:10.1016/S0309-1708(01)00026-4
.
Tolikas
P. K.
,
Sidiropoulos
E.
&
Tzimopoulos
C. D.
1984
A simple analytical solution for the Boussinesq one-dimensional groundwater flow equation
.
Water Resources Research
20
(
1
),
24
28
.
Tolikas
P. K.
,
Damaskinidou
A.
&
Sidiropoulos
E.
1990
Series representation of Flux for the Boussinesq equation
.
Water Resources Research
26
(
12
),
2881
2886
.
Tzimopoulos
C.
&
Terzides
G.
1975
Écoulement non permanent dans un sol drainé par des fosses parallèles
.
Journal of Hydrology
27
,
73
93
.
Tzimopoulos
C.
&
Tolikas
P.
1980
Technical and theoretical aspects in artificial ground water recharge
.
ICID Bulletin (International Commission on Irrigation and Drainage)
29
(
2
),
40
44
.
Tzimopoulos
C.
,
Papadopoulos
K.
,
Evangelides
C.
&
Papadopoulos
B.
2019
Fuzzy solution to the unconfined aquifer problem
.
Water
11
(
54
),
1
19
.
Tzimopoulos
C.
,
Papadopoulos
K.
,
Evangelides
C.
&
Spyrides
A.
2021
Recharging and discharging of unconfined aquifers. Case of nonlinear Boussinesq equation
. In
Proceedings of the Eighth International Conference on Environmental Management,Engineering, Planning and Economics (CEMEPE 2021) and SECOTOX Conference
,
July 20–24, 2021
,
Thessaloniki, Greece
, pp.
74
80
.
Vorobiev
D.
&
Seikkala
S.
2002
Towards the theory of fuzzy differential equations
.
Fuzzy Sets and Systems
125
,
231
237
.
doi:10.1016/S0165-0114(00)00131-7
.
Wiedeburg
O.
1890
Über die hydrodiffusion
.
Annalen der PhysiK
277
(
12
),
675
711
.
Zadeh
L. A.
1965
Fuzzy sets
.
Information and Control
8
,
338
353
.
This is an Open Access article distributed under the terms of the Creative Commons Attribution Licence (CC BY 4.0), which permits copying, adaptation and redistribution, provided the original work is properly cited (http://creativecommons.org/licenses/by/4.0/).

Supplementary data