In this paper, we propose a comprehensive methodological framework for solving the fuzzy groundwater flow problem in a simpler and faster way based on numerical analysis. In particular, a novel simplified matrix explicit inverse formula is proposed as an efficient method to solve the fuzzy Finite Difference numerical scheme, called Crank–Nicolson implicit scheme. The main advantage of the proposed method is that it offers an efficient and simpler solution to the algebraic tridiagonal system of equations that appeared in the fuzzy Crank–Nicolson scheme, without affecting the accuracy of the results. Without doubt, the simulation of an unconfined aquifer flow, using the Finite Difference Method, is often time-consuming due to the costly calculations required to solve the implicit Finite Difference schemes. This process becomes even more difficult when the physical problem is solved in its fuzzy form where the calculations become more complex and greatly increase. However, problem uncertainties are not considered negligible and must be included in the final calculations for a more rational management of the water body. This creates the need to find a new and faster way to solve the fuzzy implicit schemes, as proposed in this work. The simplified matrix explicit inverse formula was applied to the fuzzy Crank–Nicolson scheme to solve the fuzzy partial differential equation of Boussinesq and, hence, the simulation of the physical problem. Its results were compared with the corresponding results of the Thomas algorithm, which until today, is the most common method of solving the tridiagonal system of equations. The case under consideration refers to a sudden rise in the lake's water level, thus resulting in the aquifer recharging from the lake. The performance of the proposed method was more than satisfactory in terms of calculation accuracy and reliability where the numerical results are matched with the differences appearing from the 12 decimal point onward between the two examined methods. The proposed method tested also considered the running times, achieving much better times compared with the Thomas algorithm, where the differences ranged from a few seconds to several hours depending on the examined time and space steps.

• This study proposes a novel solution method of the fuzzy Crank–Nicolson scheme.

• It provides better running times for the fuzzy numerical schemes.

• It compares the proposed simplified explicit inverse formula with the Thomas algorithm.

• Faster conclusions for real physical problems contain fuzziness.

• This study provides a strong advantage to engineers for efficient water management planning.

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

### Physical problem in the fuzzy environment

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.

In the case of the aquifer illustrated in Figure 2(a), a sudden rise in the lake's water level is observed and the aquifer flow is described by the following PDE of Boussinesq:or by its linearized representation (De Ridder & Zijlstra 1994):where K is the hydraulic conductivity of the aquifer, S is the specific yield of the aquifer or drainable pore space, D is the depth of the lake, h(x, t) is the depth of the aquifer, and x, t is the coordinates (space and time).
In addition, the water balance of the aquifer could not be ignored. Although the groundwater movement takes place with very low speeds, the cross section in which it moves is huge, and, thus, the final water volumes that are transferred are very large. The following non-dimensional equations (Equations (3a) – analytical solution and (3b) – numerical solution) are used to calculate the volume of water exchanged between the aquifer and the lake (Tzimopoulos et al. 2018b):

### Solutions of the fuzzy linear PDE of Boussinesq

#### Crisp BVPs – corresponding system

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:

Figure 3

Fuzzy boundary condition .

Figure 3

Fuzzy boundary condition .

Close modal

#### Fuzzy Crank–Nicolson scheme

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.

### Solving the C–N algebraic system of equations

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.

#### TriDiagonal Matrix Algorithm (TDMA)

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 matrix is written as follows:

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

 Algorithm: Numeric algorithm for solving the tridiagonal linear system Step 1: Step 2: for ⁠, docompute and simplify   end do  Step 3: Compute the solution vector using for do end do
 Algorithm: Numeric algorithm for solving the tridiagonal linear system Step 1: Step 2: for ⁠, docompute and simplify   end do  Step 3: Compute the solution vector using for do end do

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

#### Simplified explicit matrix inverse formula

However, it is a well-known property that if a matrix is multiplied by its inverse, then the result is the unit-identity matrix ( ). Using this property, the following final relation is obtained:

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.

According to the theory of da Fonseca & Petronilho (2001), there is the classic problem of a Toeplitz matrix inverse. Assume that , and in the following matrix, with and . Then, the following tridiagonal Toeplitz matrix Σ is created:

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:

Therefore, to calculate the inverse of P, it is enough to invert the simpler matrix Z and after multiplying it by . The above can be easily proved with the following very simple mathematical proof:where I is the identity matrix.

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.

Starting from relation 6 and for , then . Finally, the formula becomes:where:and

### Numerical results of the fuzzy PDE of Boussinesq

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.

To calculate the inverse of the above matrix P, Equation (7) will be used, and then it follows that:where:After applying the SEIF in matrix Z, the following final matrix will emerge:

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

with .

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

Table 1

Numerical results from the TDMA and SEIF for the simulation time of 3 days

TDMA
x   48.875000000000000 50.000000000000000 51.125000000000000
40 46.743908912644800 47.598770605797200 48.453632298948900
80 44.805633592441600 45.414798414018700 46.023963235595600
120 43.200683191675300 43.606403596253900 44.012124000832300
160 41.988504724172200 42.240568703292700 42.492632682413100
200 41.150929618036200 41.296822104829500 41.442714591623000
240 40.618697631963400 40.697124092353100 40.775550552743000
280 40.303465834124000 40.341933334224300 40.380400834324700
320 40.120983206161600 40.136319105534200 40.151655004906900
360 40.000000000000000 40.000000000000000 40.000000000000000
SEIF
x   48.875000000000000 50.000000000000000 51.125000000000000
40 46.743908912645400 47.598770605797700 48.453632298949400
80 44.805633592442300 45.414798414019500 46.023963235596200
120 43.200683191676300 43.606403596254900 44.012124000833200
160 41.988504724173800 42.240568703294400 42.492632682414800
200 41.150929618038000 41.296822104831500 41.442714591624800
240 40.618697631965000 40.697124092354800 40.775550552744600
280 40.303465834125300 40.341933334225700 40.380400834326000
320 40.120983206162400 40.136319105535100 40.151655004907700
360 40.000000000000000 40.000000000000000 40.000000000000000
TDMA
x   48.875000000000000 50.000000000000000 51.125000000000000
40 46.743908912644800 47.598770605797200 48.453632298948900
80 44.805633592441600 45.414798414018700 46.023963235595600
120 43.200683191675300 43.606403596253900 44.012124000832300
160 41.988504724172200 42.240568703292700 42.492632682413100
200 41.150929618036200 41.296822104829500 41.442714591623000
240 40.618697631963400 40.697124092353100 40.775550552743000
280 40.303465834124000 40.341933334224300 40.380400834324700
320 40.120983206161600 40.136319105534200 40.151655004906900
360 40.000000000000000 40.000000000000000 40.000000000000000
SEIF
x   48.875000000000000 50.000000000000000 51.125000000000000
40 46.743908912645400 47.598770605797700 48.453632298949400
80 44.805633592442300 45.414798414019500 46.023963235596200
120 43.200683191676300 43.606403596254900 44.012124000833200
160 41.988504724173800 42.240568703294400 42.492632682414800
200 41.150929618038000 41.296822104831500 41.442714591624800
240 40.618697631965000 40.697124092354800 40.775550552744600
280 40.303465834125300 40.341933334225700 40.380400834326000
320 40.120983206162400 40.136319105535100 40.151655004907700
360 40.000000000000000 40.000000000000000 40.000000000000000

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.

Figure 6

Dimensional depth profiles as a function of x (m) for 10 and 120 days and for both analytical and numerical results where the SEIF is simplified.

Figure 6

Dimensional depth profiles as a function of x (m) for 10 and 120 days and for both analytical and numerical results where the SEIF is simplified.

Close modal

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.

### Aquifer's recharging water volumes

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 m3/m is ±497.52 m3/m, and this deviation could not be detected without using the fuzzy methodology.

Figure 8

Membership functions at three different times.

Figure 8

Membership functions at three different times.

Close modal

### Comparative running time results of solving the fuzzy C–N scheme

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.

Table 2

TDMA and SEIF running times (seconds) in five different simulation times and two spatial steps

Simulation time (days)Δξ = 0.01
Δξ = 0.0001
Matrix sizeRunning time (s)
Matrix sizeRunning time (s)
TDMASEIFTDMASEIF
8 × 8 0.32 0.22 800 × 800 12.98 3.14
10 23 × 23 0.46 0.25 2300 × 2300 156.05 47.85
40 43 × 43 0.85 0.49 4300 × 4300 1928.87 404.29
120 68 × 68 6.95 4.56 6800 × 6800 14,276.32 2930.10
200 98 × 98 46.51 33.55 9800 × 9800 44,020.10 9495.14
Simulation time (days)Δξ = 0.01
Δξ = 0.0001
Matrix sizeRunning time (s)
Matrix sizeRunning time (s)
TDMASEIFTDMASEIF
8 × 8 0.32 0.22 800 × 800 12.98 3.14
10 23 × 23 0.46 0.25 2300 × 2300 156.05 47.85
40 43 × 43 0.85 0.49 4300 × 4300 1928.87 404.29
120 68 × 68 6.95 4.56 6800 × 6800 14,276.32 2930.10
200 98 × 98 46.51 33.55 9800 × 9800 44,020.10 9495.14

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.

### Groundwater sustainable development and future research work

#### Sustainable Development Goals – interlinkage with groundwater

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.

#### Future work and next steps

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.

Abbaspandy
S.
&
Allahviranloo
T.
2002
Numerical solution of fuzzy differential equation
.
Mathematical and Computational Applications
7
,
41
52
.
Allahviranloo
T.
,
Gouyandeh
Z.
,
Armand
A.
&
Hasanoglu
A.
2015
On fuzzy solutions for heat equation based on generalized Hukuhara differentiability
.
Fuzzy Sets and Systems.
https://doi.org/10.1016/j.fss.2014.11.009
.
Bayrak
M. A.
&
Can
E.
2020
Numerical solution of fuzzy parabolic differential equations by a finite difference methods
.
Turkish World Mathematical Society Journal of Applied and Engineering Mathematics
10
(
4
),
886
896
.
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
(
3
),
581
599
.
https://doi.org/10.1016/j.fss.2004.08.001
.
Bede
B.
&
Stefanini
L.
2011
Solution of fuzzy differential equations with generalized differentiability using LU-parametric representation
. In
Proceedings of the 7th Conference of the European Society for Fuzzy Logic and Technology, EUSFLAT 2011 and French Days on Fuzzy Logic and Applications,
LFA 2011, pp. 785–790
.
Bede
B.
&
Stefanini
L.
2013
Generalized differentiability of fuzzy-valued functions
.
Fuzzy Sets and Systems
230
,
119
141
.
https://doi.org/10.1016/j.fss.2012.10.003
.
Boussinesq
J.
1904
Recherches theoriques sur l'ecoulement des nappes d'eau infiltrees dans le sol et sur le debit des sources
.
Journal de Mathématiques Pures et Appliquées
10
,
5
78
.
Britz
D.
,
Østerby
O.
&
Strutwolf
J.
2003
Damping of Crank–Nicolson error oscillations
.
Computational Biology and Chemistry
27
,
253
263
.
doi:10.1016/S0097-8485(02)00075-X
.
Caratelli
D.
&
Ricci
P. E
, .
2021
Inversion of tridiagonal matrices using the Dunford-Taylor's integral
.
Symmetry.
https://doi.org/10.3390/sym13050870
.
Chan
T.
1984
Stability analysis of finite difference schemes for the advection-diffusion equation
.
SIAM Journal on Numerical Analysis
21
(
2
),
272
284
.
Conte
S. D.
&
deBoor
C.
1972
Elementary Numerical Analysis
.
McGraw-Hill
,
New York
.
Crank
J.
&
Nicolson
P.
1947
A practical method for numerical evaluation of solutions of partial differential equations of the heat-conduction type
.
Mathematical Proceedings of the Cambridge Philosophical Society
.
https://doi.org/10.1017/S0305004100023197
.
De Ridder
N.
&
Zijlstra
G.
1994
Seepage and groundwater flow
. In:
Ritzema, H. P. (ed.).
Drainage Principles and Applications
, 2nd edn.
ILRI Publication 16. ILRI, Wageningen, The Netherlands, pp. 305–340
.
Encinas
A. M.
&
Jiménez
M. J.
2018
Explicit inverse of a tridiagonal (p,r)-Toeplitz matrix
.
Linear Algebra and Its Applications.
https://doi.org/10.1016/j.laa.2017.06.010
.
da Fonseca
C. M.
&
Petronilho
J.
2001
Explicit inverses of some tridiagonal matrices
.
Linear Algebra and Its Applications
325
(
1–3
),
7
21
.
https://doi.org/10.1016/S0024-3795(00)00289-5
.
da Fonseca
C. M.
&
Petronilho
J.
2005
Explicit inverse of a tridiagonal K − Toeplitz matrix
.
Numerische Mathematik
100
(
3
),
457
482
.
https://doi.org/10.1007/s00211-005-0596-3
.
Goetshel
R.
&
Voxman
W.
1986
Elementary fuzzy calculus
.
Fuzzy Sets and Systems
18
,
31
43
.
Guppy
L.
,
Uyttendaele
P.
,
Villholth
K. G.
&
Smakhtin
V
, .
2018
Groundwater and Sustainable Development Goals: Analysis of Interlinkages
.
UNU-INWEH Report Series, Issue 04
.
United Nations University
.
Hovda
S.
2016
Closed-form expression for the inverse of a class of tridiagonal matrices
.
Numerical Algebra, Control and Optimization.
https://doi.org/10.3934/naco.2016019
.
Huang
Y.
&
McColl
W. F.
1997
Analytical inversion of general tridiagonal matrices
.
Journal of Physics A: Mathematical and General.
https://doi.org/10.1088/0305-4470/30/22/026
.
Isaacson
E.
&
Keller
H. B.
1969
Analysis of numerical methods
.
The Mathematical Gazette.
https://doi.org/10.2307/3614614
.
Kasri
E. J.
,
Lahmili
A.
,
Soussi
H.
,
Jaouda
I.
&
Bentaher
M.
2021
Trend analysis of meteorological variable: rainfall and temperature
.
Civil Engineering Journal
7
(
11
),
1868
1879
.
http://dx.doi.org/10.28991/cej-2021-03091765
.
Khastan
A.
&
Nieto
J. J.
2010
A boundary value problem for second order fuzzy differential equations
.
Nonlinear Analysis, Theory, Methods and Applications
72
(
9–10
),
3583
3593
.
https://doi.org/10.1016/j.na.2009.12.038
.
Knott
G. D.
2000
Interpolating Cubic Splines
.
Birkhäuser
,
Boston, MA, USA
.
Langer
P.
2020
Groundwater mining in contemporary urban development for European spa towns
.
Journal of Human, Earth, and Future
1
(
1
),
1
9
.
http://dx.doi.org/10.28991/HEF-2020-01-01-01
.
Lee
W. T.
2011
Tridiagonal Matrices: Thomas Algorithm. MS6021, Scientific Computation
.
University of Limerick
,
Limerick, Ireland
.
Lewis
J. W.
1982
Inversion of tridiagonal matrices
.
Numerisce Mathematik
38
,
333
345
.
Mallik
R. K.
2001
The inverse of a tridiagonal matrix
.
Linear Algebra and Its Applications
.
https://doi.org/10.1016/S0024-3795(00)00262-7
.
Meurant
G.
1992
A review on the inverse of symmetric tridiagonal and block tridiagonal matrices
.
SIAM Journal on Matrix Analysis and Applications
13
,
707
728
.
Mukeru
S.
&
Mulaudzi
P. M.
2021
On the inverse problem of fractional Brownian motion and the inverse of infinite Toeplitz matrices
.
Journal of Physics Communications
5
,
1
14
.
https://doi.org/10.1088/2399-6528/ac2f91
.
Nieto
J. J.
,
Khastan
A.
&
Ivaz
K.
2009
Numerical solution of fuzzy differential equations under generalized differentiability
.
Nonlinear Analysis: Hybrid Systems
3
(
4
),
700
707
.
https://doi.org/10.1016/j.nahs.2009.06.013
.
Polubarinova-Kochina
1949
On unsteady motions of groundwater during seepage from water reservoirs
.
13
(
2
).
Pozrikidis
C.
2014
An Introduction to Grids, Graphs, and Networks
.
Oxford University Press
,
New York
.
Runkel
R. L.
&
Chapra
S. C.
1993
An efficient numerical solution of the transient storage equations for solute transport in small streams
.
Water Resources Research
.
https://doi.org/10.1029/92WR02217
.
Salih
A.
2010
Tridiagonal Matrix Algorithm
.
Department of Aerospace Engineering Indian Institute of Space Science and Technology
,
Thiruvananthapuram
.
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.
https://doi.org/10.2166/ws.2021.115
.
Shoaid
M.
,
Yang
M.
,
Liang
Y.
&
Rehman
G.
2021
Stability and deformation analysis of landslide under coupling effect of rainfall and reservoir drawdown
.
Civil Engineering Journal
7
(
07
),
1098
1111
.
http://dx.doi.org/10.28991/cej-2021-03091713
.
Smith
G. D.
1965
Numerical Solution of Partial Differential Equations: Finite Difference Method
.
Clarendon Press
,
Oxford, UK
.
Stefanini
L.
&
Bede
B.
2009
Generalized Hukuhara differentiability of interval-valued functions and interval differential equations
.
Nonlinear Analysis, Theory, Methods and Applications
71
(
3–4
),
1311
1328
.
https://doi.org/10.1016/j.na.2008.12.005
.
Stefanini
L.
,
Sorini
L.
&
Guerra
M. L.
2006
Parametric representation of fuzzy numbers and application to fuzzy calculus
.
Fuzzy Sets and Systems.
https://doi.org/10.1016/j.fss.2006.02.002
.
Tan
L. S. L.
2019
Explicit inverse of tridiagonal matrix with applications in autoregressive modelling
.
IMA Journal of Applied Mathematics.
https://doi.org/10.1093/imamat/hxz010
.
Thomas
L. H.
1949
Elliptic Problems in Linear Differential Equations Over a Network
.
Watson Sci Lab Report Columbia University
,
New York
.
Tzimopoulos
C.
,
Evangelides
C.
,
K.
&
B
, .
2016
A second order fuzzy differential equation for the case of a semi-confined aquifer
.
International Journal of New Technology and Research (IJNTR)
9
,
23
29
.
Tzimopoulos
C.
,
K.
&
Evangelides
C
, .
2018a
Fuzzy analytical solution to horizontal infiltration
. In
AIP Conference Proceedings
.
https://doi.org/10.1063/1.5043916
.
Tzimopoulos
C.
,
K.
,
Evangelides
C.
&
B.
2018b
Fuzzy solution to the unconfined aquifer problem
.
Water (Switzerland)
11
(
1
),
8
13
.
https://doi.org/10.3390/w11010054
.
Usmani
R. A.
1994
Inversion of a tridiagonal Jacobi matrix
.
Linear Algebra and Its Applications.
https://doi.org/10.1016/0024-3795(94)90414-6
.
Whorter
Mc. D.
&
D
, .
1977
Ground Water Hydrology and Hydraulics
.
Water Resources Publications
,
Littleton, CO, USA
.
Zhang
H.
,
Lu
J.
&
Hu
Q.
2014
Exponential growth of solution of a strongly nonlinear generalized Boussinesq equation
.
Computers and Mathematics with Applications
.
https://doi.org/10.1016/j.camwa.2014.10.012
.
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/).