Fuzzy numerical solution to the unconfined aquifer problem under the Boussinesq equation Open Access

The Boussinesq equation, proposed by Boussinesq (1904), is a one-dimensional nonlinear parabolic second-order partial differential equation and describes the physical problem of horizontal water flow concerning unconfined aquifers without precipitation. For this equation the following assumptions apply: (1) the inertial forces are negligible and (2) the horizontal component of velocity V_x does not vary depending on depth and its function of (x, t). Many researchers have dealt with classical (two-valued) logic and gave analytical (Boussinesq 1904; Polubarinova-Kochina 1949, 2015; Yeh 1970; Glover 1978; Koussis 1979; De Ridder & Zijlstra 1994; Zhang et al. 2014) and numerical (Neuman & Witherspoon 1971; van Schilfgaarde 1974; Tzimopoulos & Terzidis 1975; Tzimopoulos 1976; Colmenares & Neilan 2016) solutions to the equation but all of the aforementioned problems convey fuzziness (Cox 1994). In general, until recently, in practice, an attempt has been made to solve a classical mechanical problem introducing many uncertainties in the solution process.

The physical problem described by the Boussinesq equation can produce ambiguities and uncertainties regarding: the definition of the initial flow condition, the way the Boussinesq equation becomes linearized, the definition of drain spacing and the hydraulic conductivity, boundary conditions, etc. (Tzimopoulos et al. 2018a, 2018b). The uncertainties and ambiguities that groundwater physical problems present must be taken into account, and in not being included in the final calculations, can lead to incorrect conclusions and management practices, with disastrous consequences for surface and groundwater resources and therefore enormous environmental, social and economic impacts. For all the above reasons, the solution of this problem should be obtained by fuzzy algorithms.

Fuzzy logic theory is derived from the development of the fuzzy sets theory of Lofti Zadeh (Zadeh 1965) and is a powerful tool to model ambiguity, and its development gave a big boost not only to theoretical problems (Bhaskar et al. 2004; Nieto & Rodríguez-López 2006; Aminikhah 2015) but also to engineering and hydraulic problems (Guo et al. 2003a, 2003b). When we are studying in fields of physics and engineering, we often meet problems of fuzzy partial differential equations which have to be solved by numerical methods because the exact solution of these problems, as in classical logic, can be found only in some special cases (Boussinesq 1904; Polubarinova-Kochina 2015). This case of numerical methods for solving fuzzy differential equations has been rapidly growing in recent years.

During recent years, some analytical and numerical methods have been proposed in order to solve fuzzy differential equations. Initially, the concept of fuzzy derivative was introduced by Chang and Zadeh (Chang & Zadeh 1972), followed by Dubois and Prade (Dubois & Prade 1982), who used the extension principle in their approach. Fuzzy differential functions were studied by Puri and Ralescu (Puri & Ralescu 1983), who extended Hukuhara's derivative (H-derivative) (Hukuhara 1967) of a set of values appearing in fuzzy sets. Kaleva and Seikkala (Kaleva 1987, 1990; Seikkala 1987) developed the fuzzy initial value problem. But this method has 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 Bede and Gal (Bede & Gal 2005) and Stefanini and Bede (Stefanini & Bede 2009) and overcomes this drawback. This new derivative is defined for a larger class of fuzzy functions than the Hukuhara derivative. Allahviranloo et al. in Allahviranloo et al. (2015) introduced (gH-p) differentiability for partial derivatives as an extension of the above theory. Tzimopoulos et al. in Tzimopoulos et al. (2018a, 2018b) used the above method and gave a fuzzy analytical solution to a parabolic differential equation and also Tzimopoulos et al. in Tzimopoulos et al. (2018a, 2018b, 2020) gave a fuzzy analytical solution to an unconfined aquifer problem described by the Boussinesq equation.

Analytical methods approach real problems with changing parameters in space and time only in limiting cases and this disadvantage leads to the use of numerical methods such as finite difference methods transformed in the fuzzy environment, which approach any physical problem. The numerical method for solving fuzzy differential equations was introduced by Ma et al. (Ma et al. 1999). Subsequently, numerical solutions of fuzzy differential equations were examined by Friedman (Friedman et al. 1999), Bede (Bede 2006) and Abbasbandy (Abbasbandy & Allahviranloo 2002). The existence of solutions for fuzzy partial differential equations was investigated also by Buckley and Feuring (Buckley & Feuring 1999); their proposed method works only for elementary partial differential equations. Based on the Seikkala derivative, Allahviranloo (Allahviranloo 2002), and Kermani and Saburi (Kermani & Saburi 2007) use a numerical method which is an explicit difference method to solve partial differential equations. Farajzadeh (Farajzadeh 2010) gives an explicit method for solving fuzzy partial differential equations. Uthirasamy (Uthirasamy 2014) gives studies on numerical solutions of fuzzy boundary value problems and fuzzy partial differential equations. More recently, Samarinas et al. (Samarinas et al. 2018) give a fuzzy implicit scheme to solve numerically the physical problem of horizontal infiltration.

In this paper, a fuzzy numerical solution of the linearized one-dimensional Boussinesq equation, with initial and boundary conditions, is presented applying the Crank–Nicolson scheme (Crank & Nicolson 1947) in a fuzzy environment. Two different cases of the problem of unsteady flow in a semi-infinite unconfined aquifer bordering a lake are examined. In the first case there is a sudden rise of the lake's water level, thus the aquifer recharging from the lake, and for the second case a sudden drop of the lake's water level, thus the aquifer discharging to the lake. In both cases special emphasis is given to the aquifer boundary conditions, which are considered to be uncertain, as their uncertainties are the ones that affect to a greater extent the solution of the problem. Moreover, the hydraulic parameters of this problem are considered crisp as well as the geometric parameters. In general the hydraulic communication between the lake and the aquifer has an important effect on the control of the riparian ecosystem and it can also alter the water chemistry. Engineers must have knowledge of the above effects on ecological and hydrological processes for water resource management purposes.

In order to present in a comprehensive way the methodology used in this work, initially for the understanding of the physical problem and the difficulties it presents but also for the development of its solution, it is considered appropriate to present the flowchart shown in Figure 1. In this flowchart, the methodology of solving the physical problem is developed, determining that the ambiguities of the problem require the use of fuzzy partial differential equations for accurate results. Furthermore, emphasise is given to the two different approaches of solving fuzzy partial differential equations, highlighting the advantages and disadvantages of each method in order to solve the problem.

Definition 1. We denote by the class of fuzzy subsets ⁠, satisfying the following properties (Puri & Ralescu 1983; Kaleva 1987):

• 1.

  is normal, that is, there exists with ⁠;

• 2.

  is a convex fuzzy set, that is, ⁠;

• 3.

  is upper semi-continuous on ⁠;

• 4.

  is compact, where denotes the closure of A.

Then is called the space of fuzzy numbers.

Definition 2. Let be a fuzzy number. Fundamental concepts in fuzzy theory are the support, the level-sets (or level-cuts) and the core of a fuzzy number. The α-level set of (or simply α-cut) is defined by and for α = 0 it is the closure of the support (Stefanini & Bede 2014).

Definition 3. The necessary and sufficient conditions for to define a fuzzy number are as follows (Goetschel & Voxman 1986; Stefanini et al. 2006; Nieto et al. 2009):

• 1.

  is a bounded monotonic non-decreasing left-continuous function for all and right-continuous for ⁠;

• 2.

  is a bounded monotonic non-increasing left-continuous function for all and right-continuous for ⁠;

• 3.

  for α = 1, which implies ⁠.

Then it is easy to see that D is a metric in and has the following properties:

• 1.

  ⁠;

• 2.

  ⁠;

• 3.

  and is a complete metric space, for all and ⁠.

Definition 6.

A. First order

A fuzzy-valued function H of two variables is a rule that assigns to each ordered pair of real numbers (x, t), in a set D a unique fuzzy number denoted by . Let ⁠, and be real-valued functions and partially differentiable with respect to x and t. We say that (Khastan & Nieto 2010; Allahviranloo et al. 2015; Mondal & Roy 2016):

• is [(i)-p]-differentiable w.r.t x at if

  
• is [(ii)-p]-differentiable w.r.t x at if

  

Notation. The same is valid for

Definition 7.

B. Second order

Let ⁠, and be [gH-p]-differentiable at with respect to x. We say that (Khastan & Nieto 2010; Allahviranloo et al. 2015):

• 1.
  is [(i)-p]-differentiable w.r.t x if:

  
• 2.
  is [(ii)-p]-differentiable w.r.t x if:

  

The method of solving the above problem is based on the selection of the appropriate type of derivative, according to the theory of Khastan & Nieto (2010) and Bede & Stefanini (2013). With such a selection, the above problem is transformed to a corresponding system of boundary value problems.

A domain with the valid solution results, in which the derivatives have level sets according to the above theory of differentiability. Analytical solutions of the above systems can be found in Tzimopoulos et al. (2018a, 2018b).

The main idea in the finite difference method is to replace the derivative in a partial differential equation with difference quotients.

Assume that u is a function of the independent crisp variables x and t. Subdivide the x–t plane into sets of equal rectangles of sides h, k, by equally spaced grid lines parallel to ⁠, defined by and equally spaced grid lines parallel to ⁠, defined by where m and n are positive integers with and (Figure 4).

If System 2 satisfies Equation (5) then the fuzzy Crank–Nicolson scheme develops as follows.

In the cases of Systems (2,1) and (2,2) it is apparent that the fuzzy Crank–Nicolson formulas are similar to Systems (1,2) and (1,1) respectively.

Generally, when numerical methods are applied, comparison with the analytical solution of the problem – when it exists – is a way of verifying the results. In this physical problem the analytical solution exists and for this reason, the fuzzy Crank–Nicolson scheme is applied below and compared with the results of the analytical solution. Also, we should mention that the fuzzy Crank–Nicolson scheme is solved with the Thomas method (algorithm) which is a simplified form of Gaussian elimination that can be used to solve tridiagonal systems of equations (Smith 1965). Another very important point to note when using differential equations is that it is a well-posed problem, as according to Jacques Hadamard (Hadamard 1902).

We used the fuzzy Crank–Nicolson scheme (11) to approximate the fuzzy analytical solutions (10a) and (10b) for different times and points. The optimal spatial and temporal steps that emerged were dt = 0.0000018, dx = 0.01 thus r = 0.018 for the recharging case, and dt = 0.00000108, dx = 0.01 thus r = 0.0108 for the discharging case. Figures 6 and 7 illustrate the dimensional depth profiles as a function of x for the recharging and discharging cases respectively.

Figure 8 illustrates, as an example, the membership functions for the recharging case for one point in three different times. According to the membership functions, it appears that the fuzziness follows the time. This means that as the time tends to grow, higher levels of fuzziness arise.

In addition, Figures 9 and 10 illustrate in 3D mode the dimensional depth profiles as a function of x, using the fuzzy numerical solution, for the time of 120 days, presenting in detail the water movement.

As a first point of view according to the above graphs, we can observe that the fuzzy numerical results tend to approach the fuzzy analytical solution satisfactorily. However, a more reliable comparison method, which comes to confirm the visual results, is based on the transformed Haussdorff distance. This metric (Equation (2)) was applied to compare the final results and Figure 11 presents the transformed Haussdorff distances for both cases in three different times.

Based on the results, it is observed that the maximum distances between the fuzzy numbers of the two methods appeared in the results for level α-cut = 0, which is absolutely reasonable since in this α-cut level the results are extracted taking into account the full percentage of fuzziness. In addition, it is observed that the distances of the fuzzy numbers decrease as the observation times increase, which proves the stability of the fuzzy computational scheme.

The theory of fuzzy logic as well as the fuzzy partial differential equations under the generalized Hukuhara derivative proved to be valuable tools in applications for groundwater physical problems which were characterized largely by ambiguities and uncertainties. Researchers and engineers can now solve practical problems including the vital factor of fuzziness in their decision making and planning.

Analyzing the results and taking into account the values of the transformed Haussdorff distance, it was concluded that the proposed fuzzy Crank–Nicolson procedure approaches with satisfactory accuracy the fuzzy analytical results. Therefore, the proposed Crank–Nicolson scheme in a fuzzy environment is perfectly suited for scientific purposes and treatment of similar physical problems as that of this paper and also various other problems such as the enrichment of an underground aquifer by a ditch, etc.

Regarding the recharging and discharging cases, the function is [(ii)-p]-differentiable with respect to x and [(i)-p] with respect to t. The function is [(ii)-p]-differentiable with respect to x.

It must always be borne in mind that the percentage of fuzziness must be wisely defined and under no circumstances frivolous and careless as it can lead to misleading results and wrong management decisions. The expertise and experience of the engineers or researchers can contribute to the appropriate estimation of the percentage of uncertainty.

In this work, an implicit finite difference numerical scheme was proposed to solve the linearized Boussinesq equation in a fuzzy environment. An extension of this research could be the solution with both analytical and numerical approach of the non-linear fuzzy equation of Boussinesq. Also, the solution of the tridiagonal system of the fuzzy numerical scheme could be improved in terms of the computing process speed, as in a fuzzy environment more and complex calculations are required and the numerical schemes become costly in time.

No funding to declare.

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