An efficient method to solve the fuzzy Crank–Nicolson scheme with application to the groundwater flow problem

The nonlinear parabolic Partial Differential Equation (PDE) of Boussinesq proposed by Boussinesq (1904) describes the horizontal water flow of an unconfined aquifer with the following assumptions: (i) the inertial forces are negligible and (ii) the horizontal component of velocity u_x does not vary depending on depth and is a function of x and t. This problem was applied extensively to both, recharging and discharging of an unconfined aquifer, due to a sudden change in the head at origin, and to other hydraulic problems. The solution of the linearized form of the Boussinesq equation using the Finite Difference (FD) method, such as the Crank–Nicolson (C–N) implicit scheme (Crank & Nicolson 1947), is obtained by transforming the above PDE linearized equation in a system of the simultaneous linear algebraic equations with non-zero coefficients only on the main diagonal, the lower diagonal, and the upper diagonal, and this is called the tridiagonal system of equations (Salih 2010). Such schemes provide unconditional stability (Isaacson & Keller 1969; Chan 1984), producing reliable and accurate results, but their solution is quite complex and requires the use of appropriate algorithms and computers.

Until today, the most common method to solve a linear tridiagonal system of equations is the Thomas algorithm, also called the TriDiagonal Matrix Algorithm (TDMA) (Thomas 1949), which is essentially the result of applying the Gaussian elimination to the tridiagonal system, and it works efficiently with these systems. Although the TDMA is the most simple and economical in terms of memory requirements and the computational time method (Bieniasz 2001), sometimes this method presents certain drawbacks in large-scale problems involving a large number of variables (Runkel & Chapra 1993). It can also be unstable, in some rare cases, under a specific condition for a singular or non-singular tridiagonal matrix (Lee 2011). Mainly, difficulties are found in storing the coefficients and managing the large-scale matrices, resulting from the discretization process. A great amount of memory storage is required to solve such a system in combination with the already large number of matrices used by the computational scheme, thus creating difficulty and, in many cases, prohibitive conditions for the whole process. Particularly, if the C–N scheme is to be used in its fuzzy form, the computational needs for its solution should be greatly increased. Thus, it is fully acceptable that the TDMA has many advantages against other methods, when dealing with physical problems in a classical-binary logic, but lags in complex and large-scale computational processes of the fuzzy environment. This fact has created the need to find a new and faster way to simplify and improve the solution of the tridiagonal system, so that the fuzzy FD scheme will be easily applicable and clearly faster.

Even though the TDMA is prevailing in solving tridiagonal systems, these systems can be solved alternatively using the property of matrix inverse. By inverting the proper tridiagonal matrix, which appears in the scheme, the system can be solved with a relatively simple and fast way. However, this does not mean that the classic methods of matrix inversion (such as ready functions in MATLAB and MATHEMATICA) can achieve better computational speeds than the TDMA. However, using the explicit inverse methods, significant changes occur due to the simplicity of the calculations that these methods present.

In addition to the FD method, the tridiagonal linear systems often occur in solving interpolation problems (Knott 2000) and boundary value problems (BVPs) (Pozrikidis 2014). For this reason, in recent years, the inverse of tridiagonal matrices has been extensively studied in different research fields, such as mathematics, engineering, physics, and more. Because of this, many numerical methods have been developed to give inverse types of such matrices (Lewis 1982; Mallik 2001; Hovda 2016), where some prove to be less and others more effective. However, explicit inverses are known only in some cases, where the tridiagonal matrix is symmetric with constant diagonals and some constraints. A detailed review is given by Meurant (1992), who also presents an explicit inverse way of a Toeplitz tridiagonal matrix by analytically solving the iterations of the Chebyshev polynomials. Regarding the inverse of some general tridiagonal matrix, different approaches have been proposed by several authors (Usmani 1994; Huang & McColl 1997). Using and extending Meurant's theory, da Fonseca & Petronilho (2001, 2005) provide an explicit inverse method of the k-Toeplitz matrix, thereby generalizing some already known inversion methods. Also, an explicit inverse method of a (p, r)-Toeplitz tridiagonal matrix is given by Encinas & Jiménez (2018) using the Schrödinger operator and the Chebyshev functions. Furthermore, Tan (2019) presents an explicit inverse way of a symmetric tridiagonal matrix, but with the constraint that the first and last elements of the diagonal should be different from the rest. Mukeru & Mulaudzi (2021) obtained an exact formula that expresses the elements of the inverse in terms of the coefficients of the inverse Szegö function associated with a Toeplitz matrix. More recently, Caratelli & Ricci (2021) present an expression for the inverse of a general non-singular complex-valued tridiagonal matrix, based on Dunford–Taylor's integral, giving examples of the tridiagonal Toeplitz matrix inversion of third, fourth, and fifth orders. At this point, it should be emphasized that although matrix explicit inversion formulas are presented in the literature, their tests are often limited to simple examples of an easy matrix inversion. The application of explicit formulas to large-scale real physical problems with many variables and big simulation times, as well as a comparison with different methods of solving tridiagonal systems, is something that has not been done and is worth considering, as is shown in the present work.

Regarding the groundwater flow, the correct simulation of the problem, based on the Boussinesq equation, requires a knowledge of key parameters such as the hydraulic conductivity, specific yield, and the boundary conditions. However, it is known that an unconfined aquifer system is a complex and open system that is affected by meteorological, soil, and geological conditions as well as by human activities, and has been recently studied by researchers (Langer 2020; Kasri et al. 2021; Shoaid et al. 2021). Particularly, for the definition of boundary conditions, the exact values are often difficult to obtain mainly due to the difficult location of groundwater systems, the lack of field data observations, and the large area they cover. For this reason, ambiguities and uncertainties are created in the definition of the boundary conditions and consequently of the whole problem. These uncertainties can arise in a variety of ways such as measurement and instrument errors, the measurement methods that were used, errors in raw data processing, as well as many other errors in the process of determining aquifer's boundary conditions due mainly to human mistakes that lead to an overestimation or underestimation of the real physical condition. Because of this, and as the sensitivities in the results are significant, we pick to include, in the final calculations, the uncertainties that appear in the boundary conditions that are important to be considered to enable engineers to arrive at correct conclusions and apply optimal management practices.

In many cases, finding the solution of the described PDE of Boussinesq through numerical methods is the only effective way because the disadvantage of analytical solutions to approach real problems with changing parameters in space and time creates a gap in finding exact solutions every time, except in some special cases (Polubarinova-Kochina 1949; Zhang et al. 2014). For the fuzzy PDE of Boussinesq, the appropriate approach is to use the generalized Hukuhara derivative (Bede & Gal 2005; Stefanini & Bede 2009) and the extension of this theory regarding the partial derivatives (Allahviranloo et al. 2015;Tzimopoulos et al. 2016, 2018b). As in the case of classical logic, the same is valid also for the fuzzy PDE of Boussinesq, where the use of either fuzzy analytical solutions (Tzimopoulos et al. 2016, 2018b) or fuzzy numerical schemes (Abbaspandy & Allahviranloo 2002; Samarinas et al. 2018, 2021; Bayrak & Can 2020) to approach the problem is required.

In this paper, we are motivated by the explicit inverse formula proposed by da Fonseca & Petronilho (2001) and we are simplifying it in a way that makes it more efficient and simple when applied to the demanding fuzzy C–N scheme without sacrificing accuracy. The fuzzy C–N scheme is proposed by our previous work (Samarinas et al. 2021), where the solution of the scheme was found by using the TDMA. However, it was observed that the process was time-consuming, and an extension of the previous research is the present work applying a novel solution method to the scheme, which is validated with the results of our previous work under the same groundwater physical problem.

In short, the major contributions of this work are as follows:

• It explores the use of a proposed novel simplified explicit matrix inverse formula as an efficient method for solving the algebraic tridiagonal linear system of equations that appears in the fuzzy C–N FD implicit scheme. The simplified formula comes after modifying the tridiagonal matrix, of the implicit scheme, that is imported in the algorithm for inversion.

• It compares the proposed method with the TDMA, the traditionally used approach, under a real physical problem.

• It explores the use of the fuzzy C–N scheme, which is solved by the simplified explicit inverse methodology, to solve the unconfined aquifer's real physical problem. The case under consideration concerns the hydraulic communication of an aquifer bordering a lake and specifically the case of the aquifer's recharging from the lake due to the sudden rise in the lake's water level.

• It provides accurate, reliable, and faster results to the fuzzy problem and, thus, paves the way for the groundwater problem simulation while supporting the 2030 Sustainable Development Agenda.

The methodological approach applied in the present work is presented in the following flowchart (Figure 1) for the convenience of the reader in an explanatory form. In essence, it focuses on five main steps to achieve the final goal, which is the fastest solution of the fuzzy PDE of Boussinesq with the FD method. Using an efficient method of solving the C–N algebraic tridiagonal system of equations, based on the proposed simplified explicit matrix inverse formula, was the crucial step for obtaining a faster solution of the overall fuzzy physical problem of groundwater flow.

The physical problem of the horizontal water flow concerns unconfined aquifers without precipitation. It is described by the one-dimensional nonlinear parabolic second-order PDE, known as the Boussinesq equation. The physical problem can produce ambiguities and uncertainties regarding a numerous cause, such as the definition of the initial flow condition, the way of the Boussinesq's equation becomes linearized (Tzimopoulos et al. 2018a, 2018b) or incorrect use of the other hydraulic parameters of the problem. The examined fuzzy physical problem has already been described extensively, and here, the method of transforming the PDE of Boussinesq into a fuzzy environment is briefly presented (Equation (1)). The reader is referred to Tzimopoulos et al. (2018b) for further details.

Analytical solutions of the fuzzy Equation (1) with the boundary and initial conditions (2) can be found in Tzimopoulos et al. (2018b) utilizing the theory of Khastan & Nieto (2010), Bede & Stefanini (2011, 2013), and Allahviranloo et al. (2015), by translating the above fuzzy problem to a system of second order of crisp BVPs, which are called a corresponding system of the fuzzy problem (Khastan & Nieto 2010).

In our case, the appropriate system of the fuzzy problem is as follows:

where ⁠.

Figure 4 presents the finite computational grid of the C–N implicit scheme. Based on the grid, we can observe that the difference approximations of the C–N scheme are developed at the midpoint of the time increment. The derivatives in space can be determined at the midpoint by averaging the difference approximations at the beginning and at the end of the time increment.

where ⁠.

More descriptive and are column matrices that contain the approximate solutions to be found, and the matrix P is strictly a tridiagonal matrix that contains the unknown coefficients, Q is the matrix that contains the known values of the function at the previous time steps and has a strictly tridiagonal structure as well, and are column matrices that contain the known values, and is also a column matrix that contains the condition values.

If the discretization problem consists of N − 1 internal nodes in terms of spatial dimension, matrices and Q will be (N − 1) × 1, while matrix P will be (N − 1) × (N − 1). The system would be complicated to solve, due to the number of nodes required for sufficient method accuracy, but it is observed that matrix P is sparse, i.e., it contains relatively few non-zero elements, because each node depends only on neighbor nodes. Particularly, for differential equations of the second order, such as the examined Boussinesq equation, it is proved that it is tridiagonal, i.e., its non-zero elements are gathered only in the main diagonal, the upper diagonal, and the lower diagonal. This results in less memory being stored as only 3 × (N − 1) elements are stored instead of (N − 1) × (N − 1) as defined in matrix P.

The tridiagonal matrices used in the C–N equations are usually strictly diagonally structured, the determinant of these matrices is not zero, and the matrix P is revertible. Thus, the system is not impossible, nor does it have an infinite number of solutions.

For comparing the two different examined methods of solving the fuzzy C–N scheme, it is considered necessary by the authors to present them with a detailed description, initially, for the convenience of the reader, but also for a more reliable evaluation of the methods in terms of solving the scheme and comparison in general.

As already mentioned, the implicit schemes are solved as a linear algebraic tridiagonal system of equations. The most common method to solve those systems is with the TDMA which is based on LU decomposition.

The transformed system is easy to solve by backward substitution. The first linear system can now be solved using the following algorithm:

As it is easily understood, the above algorithm consists of a forward sweep and a backward sweep, making it a double sweep algorithm.

Therefore, the linear algebraic tridiagonal system of equations resulting from the implicit C–N scheme can be solved based on the above relation, i.e., the inverse of the matrix P must be calculated.

The above inverse formula is very simple and clear and more easily to apply than the TDMA, but in order to further improve it and to achieve a better performance in the fuzzy C–N scheme, we modified the matrix P that is inverted, and we worked in the following way:

Consequently, the initial matrix A is simplified further and it now has a lower and upper diagonal −1 but also a strong main diagonal. In this way, formula 6 is further simplified.

The fuzzy C–N scheme (Equations (4a) and (4b)) is used in order to solve the fuzzy PDE of Boussinesq (Equation (1)). The linear tridiagonal system of equations is solved using the TDMA and also using the simplified explicit inverse formula (SEIF). Both methods were tested in a simulation of 3, 10, 40, 120, and 200 days and with two different steplengths (Δξ = 0.01 and 0.0001).

The SEIF is proposed in this work, while the TDMA has already been successfully tested when compared with the analytical solution of the fuzzy physical problem (Samarinas et al. 2021). It is now considered a more objective proposition if the results of the TDMA should be taken as the reference results on which to compare the corresponding results of the proposed method, in terms of accuracy and reliability.

The created matrix P of the tridiagonal system of equations is of vital importance where, in order to solve the system, this matrix must be inverted. For this reason, as an example, we will examine the case of 3 days with steplength Δξ = 0.01, time step Δτ = 0.0000018, and for the a-cut = 0. The time of 3 days is relatively short regarding the progress of the aquifer flow but was chosen as an example. Indeed, if we consider that the simulation time increases while the steplength decreases, the tridiagonal system of equations becomes much more complex with larger matrix sizes and more information in terms of the aquifer flow. In the context of these pages, it is impossible to present them in full numerical form. However, for the purpose of comparing the two methods, it is necessary to present numerical results, and the 3-day simulation fits this reason perfectly.

The above matrix was used to solve the tridiagonal system of equations in the fuzzy C–N scheme and the system is also solved with the TDMA. Therefore, a comparison of the results is now possible, and for convenience, they are defined by the following fuzzy triangular numbers (Goetshel & Voxman 1986; Stefanini et al. 2006; Nieto et al. 2009):

• the TDMA results

• the results from the SEIF

⁠

Table 1 presents the numerical results of the two examined methods in an analytical form and with 15 digits of accuracy.

As a first point of view, from the data of the above table, it is easily understood that the numerical results of the SEIF approach with very good accuracy the corresponding TDMA results. Specifically, the accuracy breaks to 12-decimal points. There is also a complete identification of the results in the boundary points 0 and 360 m, which is completely reasonable, as in those points are the known values of the initial and boundary conditions of the aquifer. In addition, the following graph (Figure 5), in which the absolute differences of the results have been taken into account, confirms to a large extent the numerical analytical results of the above table. It is easily obtained that approximately in the middle length points (160, 200, and 240 m), we have the largest differences.

Figure 6 illustrates the dimensional depth profiles for 10 and 120 days as a function of x (m) for the physical problem under consideration using both numerical (Equations (4a) and (4b)) and analytical results. The numerical approach of all solutions converges to a great extent, and for this reason, detailed differences cannot be distinguished in the context of a graphical representation. However, an illustration of the aquifer flow should be presented to draw conclusions. From the physical point of view and based on the below graph, the depth profile for 10 days approaches 900 m of the aquifer's length, while for 120 days, the depth profile goes even further, reaching 2800 m.

Considering the above, we can state that the results of the SEIF are identical to those of the TDMA, which proves the accuracy and reliability of the proposed solution approach. Moreover, the perfect accuracy of the results in all positions of the physical problem further enhances its application with longer simulation times. The same procedure was followed for the rest of the simulations of 10, 40, 120, and 200 days, where a very good accuracy was also achieved.

From a hydraulic point of view, Figure 7 shows the recharging water volumes that enter the aquifer from the lake using analytical (Equation (3a)) and numerical (Equation (3b)) solutions and both the TDMA and the SEIF. It is obvious that the results coincide to a great extent for both analytical and numerical methods and also between the TDMA and the proposed SEIF.

Furthermore, the fuzzy methodology offers a big advantage where using the membership functions the variations in the results from the crisp value can be easily and quickly identified. In this regard, Figure 8 presents the membership functions of water volumes at three different times. Observing the membership functions, it is obvious that the ranges of deviations are significant. For example, the time of 120 days, the water volume deviation of the crisp value 3316.79 m³/m is ±497.52 m³/m, and this deviation could not be detected without using the fuzzy methodology.

So far, the proposed SEIF has been tested in terms of accuracy and reliability achieving excellent results. However, the main advantage of the method is focused on the computational process speed due to the simplicity of the calculations it presents. At this point, it should be emphasized that when the FD implicit schemes are applied, tests with different space and time steps are inevitable until it achieves the desired results. This process is especially for large-scale physical problems, which are also examined in the fuzzy environment, and is time-consuming, often creating prohibitive conditions for the algorithm's applications. Therefore, approximate or even incorrect solutions with space and time steps that are approximately accurate are accepted to describe the physical problem under consideration.

For the above reasons, the proposed method was considered appropriate to test by also considering the running time in order to solve the fuzzy PDE of Boussinesq using the fuzzy C–N scheme. For the purpose of comparison, the simulation times of 3, 10, 40, 120, and 200 days with two different stepslengths (Δξ = 0.01, 0.001) were taken into account, and Table 2 presents the running times of the TDMA and also the corresponding ones of the SEIF. The reason behind the choice of small space and time steps and consequently the more detailed discretization grid of the aquifer initially focuses on the oscillatory response, that the C–N implicit scheme has, to sharp initial transients in the real physical processes of the aquifer. The damping of these oscillations can be avoided using relatively small values of the r coefficient (r = 0.018 in our case) in the C–N scheme (Britz et al. 2003). Additionally, the need to compare the two examined methods in terms of running times requires small steps so that the algorithms are more computationally expensive and the comparison as objective as possible. Of course, the two methods were tested on the same computer with processing power Intel(R) Core (TM) i7-4820 K CPU @ 3.70 GHz and 19 GB of RAM memory using MATLAB R2017b to build and run the algorithms.

It is obvious that the running times for the lengthstep Δξ = 0.01 are relatively close. However, when the lengthstep gets smaller and goes to Δξ = 0.0001, the running times have significant differences. The ascendancy of the SEIF for solving the tridiagonal system of equations and consequently the simulation of the physical problem at different times are unquestionable, thus creating a fertile ground for the large-scale simulation of physical problems in a fuzzy environment using the FD methods.

In general, water is strictly addressed by the United Nations 2030 Development Agenda and its Sustainable Development Goals (SDGs). However, the SDGs of the Agenda do not, as a rule, account explicitly for the significant role that groundwater plays and will continue to play in sustainable development. Undoubtedly, groundwater has a fundamental role to play in efforts to ensure irrigation needs and drinking water supply. It, therefore, supports food security and energy production to a great extent.

The findings of this work can contribute to the sustainable development of groundwater as the proposed methodology can provide reliable, accurate and fast results on the groundwater flow problem, including its uncertainties. In this way, it can contribute to decision-making and planning while supporting the SDGs 6 (clean water and sanitation), 12 (responsible consumption and production), and 13 (climate action), which, according to a recent report (Guppy et al. 2018), have the closest links with groundwater. Furthermore, the calculation of the groundwater recharging volumes as proposed in this work can contribute to SDG targets such as 2.4 (sustainable agriculture), 6.6 (ecosystems and water), and 15.1 (terrestrial ecosystems), which are based on groundwater availability and good management.

Further work should focus on the solution of the fuzzy C–N scheme with other solving methods that could be based on explicit inverse formulas or not and their subsequent comparison in terms of accuracy and running time with application to the real physical problem. Another interesting effort may be to develop a methodology on how to also include the uncertainties from the other Boussinesq equation parameters such as hydraulic conductivity in the final calculations. In addition, further research could be to solve the linearized Boussinesq equation with other numerical approaches such as the finite-element method considering the problem of uncertainties. Furthermore, the solution of the nonlinear fuzzy Boussinesq equation with both analytical and numerical methods would be scientifically important. Last but not least, it should be mentioned that in the present work, we treated an ideal aquifer with an assumption that the area has homogenous soil. Available observation point data can be used in future research for validation needs.

In this work, we present the simplified explicit matrix inverse formula as an efficient method for solving the fuzzy C–N implicit scheme. In particular, we focus on the development of a methodological framework that includes the ambiguities and uncertainties of the unconfined aquifer's physical problem, in the final calculations, in a reliable and faster way than has been done so far with other approaches. As demonstrated, the proposed method produced noticeable gains in accuracy and computational cost-efficiency in all the example simulation times compared with the TDMA traditional approach with a satisfactory level of accuracy and its better performance in terms of calculation speed.

In the present case, solving the fuzzy C–N scheme with the TDMA, we noted that for long simulation times, coupled with small lengthsteps, the process was time-consuming even for very powerful computers. This fact made its use difficult, but in no case could it be a disadvantage or even a reason for not using it. By applying the proposed SEIF, we show that the fuzzy implicit scheme can now be easily solved for delivering important information of the real physical problem. The inclusion of the ambiguities and uncertainties of the unconfined aquifer physical problem is crucial in terms of conclusions about the aquifer's flow and leading to the right decisions for its sustainable management. Our vision is that this methodology can support researchers and engineers to have faster conclusions for the fuzzy physical problem and taking more reliable decisions for aquifer's management and planning.

The groundwater flow physical problem presents a series of ambiguities, but we consider that the ambiguities in the boundary conditions are the ones that affect the results to a greater extent and, therefore, in this work, an ambiguity of 15% was chosen in the boundary conditions of the aquifer. However, in order to provide reliable results from this methodology, some limitations must be overcome. In general, in problems that convey fuzziness, the contribution of the researcher is considered important as the degree of uncertainty that will be given to the respective parameters should be representative of the problem. This is mainly achieved by the experience and expertise of the researcher or engineer with the real physical problem. Furthermore, it should be highlighted that the proposed simplified explicit matrix inverse formula can work only with tridiagonal Toeplitz matrices. Therefore, it cannot be applied to physical problems described by system types different from the tridiagonal systems of equations.

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