Adnan Altay Altınörs
Middle East Technical University, Ankara, Türkiye.
DOI: 10.4236/apm.2026.167025   PDF    HTML   XML   5 Downloads   27 Views   Citations

Abstract

Ensuring stability and accuracy is essential in numerical modeling of subsurface flow in fractured and heterogeneous media. This study investigates one-dimensional transient groundwater flow in a finite fractured confined aquifer using a double-porosity continuum approach, for which a finite difference scheme is developed employing the Crank-Nicolson method. The validity of the proposed scheme is assessed through detailed comparisons with an established analytical solution, demonstrating close agreement and confirming the accuracy of the numerical solution. In addition, von Neumann stability analysis is employed within a systematic framework to evaluate the stability of the scheme and to provide insight into its behavior across governing parameters. The stability condition is satisfied. The results show that the analytical and numerical solutions of the finite difference equations for the fractures exhibit the closest agreement for a specific value of $\alpha =\frac{\Delta \theta }{\Delta {\lambda }^2}$ for all cases considered, where the amplification factor, G, attains its minimum.

Keywords

Validation, Accuracy, Numerical Solution, Non-Darcian Flow, Fractured Aquifer

Share and Cite:

1. Introduction

References [1] and [2] investigated one-dimensional transient groundwater flow in a finite fractured confined aquifer using a double-porosity conceptual model based on a continuum approach, and presented a numerical solution. [3] derived an analytical solution for the same problem.

A general description of the problem is presented in Figure 1. A stream forms one boundary of the aquifer, while an impervious layer forms the other. When the stream level rises sharply, it increases hydraulic pressure or piezometric head (buildup) in the aquifer, causing unsteady flow conditions.

Ensuring the stability and accuracy of numerical solutions is a fundamental requirement in the analysis of subsurface flow, particularly in fractured and heterogeneous media. The main purpose of this study is to validate the numerical solution for the problem through detailed comparisons with an established analytical solution, and adopting a systematic and quantitative framework, wherein von Neumann Stability Analysis and Discrete Perturbation Analysis are employed to evaluate the stability characteristics of the scheme and to provide deeper insight into the observed numerical behavior across the range of governing parameters. This combined approach facilitates a rigorous assessment of the reliability and accuracy of the numerical solution.

2. Theoretical Background

When dealing with the flow in a fractured porous medium, the microscopic flow patterns inside the individual pores or fractures are overlooked and it is considered that some fictitious average flow is taking place. Therefore, the concept of continuum is employed [4].

In the double-porosity conceptual model, the fractured medium is represented by two fully overlapping continua; one corresponding to the porous matrix (low-permeability blocks with primary porosity) and the other representing the fracture network (high-permeability fractures with relatively small storage volume). The fractures provide high conductivity conduits amenable to rapid hydraulic flows, whereas the high porosity blocks contain the majority of the storage [5]. Both in the fractures and in the blocks, hydraulic flow takes place. There is also an exchange of fluid between the fractures and the blocks.

Figure 1. Idealized flow system [3].

At each geometric point of a fractured porous medium, following the continuum concept and according to the double-porosity conceptual model, two sets of parameters should be introduced, the first being for the fracture flow and the second for the flow in blocks [3] [6]. In this respect, two sets of governing equations are required to describe the flow in a fractured porous medium, one for each type of porosity. These equations are then linked by a fluid exchange, or interaction, term describing the quasi-steady state transfer of fluid from fractures to blocks or vice versa in terms of difference in piezometric heads between the blocks and the fractures.

Conventional single-porosity flow models generally fail to adequately describe flow behavior in fractured aquifers. Therefore, flow within the fractures is described using the non-linear Forcheimer’s equation, since the high velocities in fractures can lead to non-Darcian flow conditions. In contrast, flow within the porous matrix blocks is assumed to follow Darcy’s law. The exchange of water between fractures and matrix blocks is modeled as a quasi-steady-state transfer process between the two continua. The aquifer characteristics and the nonlinearity determine the flow behavior [1].

2.1. Governing Equations

The governing equations for flow in both the fractures and the matrix blocks are obtained by combining the equation of motion with the continuity equation. For the matrix blocks in a fractured confined aquifer, the resulting differential equation, derived from the continuity equation together with Darcy’s law, is given by:

$B{K}_b\frac{{\partial }^2{h}_b}{\partial x^2}=S_b\frac{\partial h_b}{\partial t}-{\nu }_d$ (1)

where h_b is the piezometric head in blocks; t and x are time and space variables respectively; B is the thickness of aquifer; K_b is the hydraulic conductivity of blocks; S_b is the storage coefficient of blocks; and ν_d is the fluid transfer rate. The procedure for the derivation of this equation may be found in [7].

For the non-Darcian fracture flow, Equation (1) should be modified by replacing Darcy’s law with the Forcheimer’s equation [2].

The Forcheimer’s equation, which can be characterized as the empirical modification of Darcy’s model [7], can be expressed in terms of hydraulic gradient in the vectorial form as;

$\nabla h=-F\left(q\right)q$ (2)

where F(q) is a scalar function of the magnitude of the specific discharge vector $q$ at any point. For one-dimensional flow in the x-direction, Forcheimer’s equation is

$F\left(q\right)=aq+b{q}^2$ (3)

$-\frac{\text{d}h}{\text{d}x}=F\left(q\right)q$ (4)

where q is the magnitude of specific discharge in x-direction; a and b are the Forcheimer parameters.

Accordingly, the differential equation governing flow in fractures can be written in the following form [2]:

$\frac{B}{F\left(q\right)}\frac{{\partial }^2{h}_f}{\partial x^2}=S_f\frac{\partial h_f}{\partial t}+{\nu }_d$ (5)

where h_f is the piezometric head in fractures and S_f is the storage coefficient of fractures.

2.2. Fluid Transfer between Fractures and Blocks

Reference [3] defined ν_d as the volume of water transfer per unit horizontal area per unit time. This parameter depends on the piezometric head difference between the blocks and fractures, as well as on the geometric and hydraulic characteristics of the fractured medium.

${\nu }_d=\epsilon T_b\left(h_f-h_b\right)$ (6)

where

$\xi =\epsilon T_b$ (7)

ξ is the interaction term which incorporates the geometry and other properties of the aquifer; ε is the fluid transfer parameter; T_b is the transmissivity of the blocks defined as:

$T_b=B{K}_b$ (8)

It should be noted that the fluid transfer rate between blocks and fractures, ν_d, serves as the coupling term between the governing Equation (1) and Equation (5).

2.3. Assumptions

Non-Darcian fracture-to-matrix flow may occur and, in general, should be taken into consideration. However, [8] concluded that the impact of the non-Darcian fracture-to-matrix flow can be ignored for almost all practical purposes in a typical fractured aquifer.

Equation (1) and Equation (5) are based on the following assumptions:

1) Darcy’s law is valid for the flow in the blocks.

2) The flow in fractures is non-Darcian and governed by Forcheimer’s equation.

3) Fractures and blocks are homogenous and isotropic.

4) Thickness of the aquifer is constant.

5) Length of the aquifer is constant.

6) Flow occurs only in x-direction.

7) The geometry of fractures is unaffected by chemical dissolution or deposition.

8) The aquifer is confined and non-leaky.

9) The flow is fully saturated.

In fractured aquifers, generally the conducting capability of the blocks is much less important than their storage capability. Therefore, as an assumption, the term on the left-hand side of Equation (1) can be neglected [9]-[11]. The analytical solution of [3] is based on this assumption. However, in that solution, the flow in both the fractures and the blocks is assumed to obey Darcy’s law.

2.4. Initial and Boundary Conditions

It is assumed that the water levels in both the stream and the aquifer are initially horizontal and that no flow exists within the aquifer at the initial state. Considering a step rise in the stream stage, the corresponding initial and boundary conditions governing the piezometric head in both the fractures and the matrix blocks are expressed as follows:

Initial conditions;

$h_b=h_{0}0\le x\le Lt=0$ (9)

$h_f=h_{0}0\le x\le Lt=0$ (10)

Boundary conditions at the stream side;

$h_b=h_0+s_{0}x=0t>0$ (11)

$h_f=h_0+s_{0}x=0t>0$ (12)

Boundary conditions at the impervious side;

$\partial h_b/\partial x=0x=Lt>0$ (13)

$\partial h_f/\partial x=0x=Lt>0$ (14)

2.5. Non-Dimensionalization of the Governing Equations and Initial and Boundary Conditions

To simplify the mathematical formulation and generalize the results independently of specific units and scales, the governing equations together with the initial and boundary conditions are expressed in dimensionless form using the following variables, similar to those introduced by [3] [10], may be defined as:

$z_b=\frac{h_b-h_0}{s_0}$ (15)

$z_f=\frac{h_f-h_0}{s_0}$ (16)

where z_b and z_f are dimensionless build-ups for blocks and fractures respectively.

$\lambda =\frac{x}{L}$ (17)

where λ is the dimensionless space variable, and

$\theta =\frac{T_{f}t}{S_f{L}^2}$ (18)

where θ is the dimensionless time variable.

Equation (19) and Equation (20) represent the non-dimensional forms of the governing Equation (1) and Equation (5), corresponding to flow in the blocks and the fractures respectively.

$\frac{{\partial }^2{z}_b}{\partial {\lambda }^2}=\frac{S_b}{S_f}\frac{T_f}{T_b}\frac{\partial z_b}{\partial \theta }-L^2\epsilon \left(z_f-z_b\right)$ (19)

$\frac{{\partial }^2{z}_f}{\partial {\lambda }^2}=\frac{F\left(q\right)}{B}{T}_f\frac{\partial z_f}{\partial \theta }+\frac{F\left(q\right)}{B}{L}^2\epsilon T_b\left(z_f-z_b\right)$ (20)

where L is length of the aquifer and T_f is transmissivity of fractures.

The non-dimensional forms of the initial and boundary conditions are given as follows:

Initial conditions;

$z_b\left(\lambda ,0\right)=0$ (21)

$z_f\left(\lambda ,0\right)=0$ (22)

Boundary conditions at the stream side;

$z_b\left(0,\theta \right)=1$ (23)

$z_f\left(0,\theta \right)=1$ (24)

Boundary conditions at the impervious side;

$\frac{s_0\partial z_b}{L\partial \lambda }=0\frac{\partial z_b}{\partial y}=0$ (25)

$\frac{s_0\partial z_f}{L\partial \lambda }=0\frac{\partial z_f}{\partial \lambda }=0$ (26)

3. Methodology

3.1. Numerical Solution

The dimensionless governing Equation (19) and Equation (20) correspond to the one-dimensional parabolic model equation, i.e. the simple diffusion equation, except for the coupling term.

The model equation can be solved numerically by the Crank-Nicolson implicit scheme, which is simply the trapezoidal method adapted to the context of parabolic partial differential equations [12]. This method is of order two in space, implicit in time, unconditionally stable and has higher order of accuracy [13].

The Crank-Nicolson scheme is always second order accurate in both Δt and Δx: This means that the extra stability properties of the scheme can be exploited to take larger time steps, with for example Δx = O(Δt), and because then the truncation error is O(Δt)², good accuracy can be achieved economically [14].

The discrete Crank-Nicolson system is obtained by applying central differencing in space and averaging the temporal terms between two consecutive time levels, tⁿ and tⁿ⁺¹. The resulting algebraic equations form a tridiagonal linear system. A time-marching algorithm proceeds sequentially from the initial condition by computing the dependent variable at the new time level tⁿ⁺¹ using the known solution at the previous time level tⁿ. This process is repeated until the final simulation time is reached.

Using the Crank-Nicolson method, Equation (27) and Equation (28) are obtained as the discretized forms of the dimensionless governing Equation (19) and Equation (20) [2].

$\begin{array}{l}-\left(\frac{S_f}{2S_b}\frac{T_b}{T_f}\frac{\Delta \theta }{\Delta y^2}\right)z_b{}_{i+1}^{n+1}+\left(1+\frac{S_f}{S_b}\frac{T_b}{T_f}\frac{\Delta \theta }{\Delta y^2}\right)z_b{}_i^{n+1}-\left(\frac{S_f}{2S_b}\frac{T_b}{T_f}\frac{\Delta \theta }{\Delta y^2}\right)z_b{}_{i-1}^{n+1}\\ =z_b{}_i^n+\frac{S_f}{2S_b}\frac{T_b}{T_f}\frac{\Delta \theta }{\Delta y^2}\left(z_b{}_{i+1}^n-2z_b{}_i^n+z_b{}_{i-1}^n\right)+\Delta \theta \frac{S_f}{S_b}\frac{T_b}{T_f}{L}^2\epsilon \left(z_f{}_i^n-z_b{}_i^n\right)\end{array}$ (27)

$\begin{array}{l}-\left(\frac{B}{2F{\left(q\right)}_i^n}\frac{1}{T_f}\frac{\Delta \theta }{\Delta y^2}\right)z_f{}_{i+1}^{n+1}+\left(1+\frac{B}{F{\left(q\right)}_i^n}\frac{1}{T_f}\frac{\Delta \theta }{\Delta y^2}\right)z_f{}_i^{n+1}-\left(\frac{B}{2F{\left(q\right)}_i^n}\frac{1}{T_f}\frac{\Delta \theta }{\Delta y^2}\right)z_f{}_{i-1}^{n+1}\\ =z_f{}_i^n+\frac{B}{2F{\left(q\right)}_i^n}\frac{1}{T_f}\frac{\Delta \theta }{\Delta y^2}\left(z_f{}_{i+1}^n-2z_f{}_i^n+z_f{}_{i-1}^n\right)-\frac{L^2\epsilon T_b}{T_f}\Delta \theta \left(z_f{}_i^n-z_b{}_i^n\right)\end{array}$ (28)

Equation (5) contains the non-linear term, $\frac{B}{F\left(q\right)}\frac{{\partial }^2{h}_f}{\partial x^2}$ , where F(q) depends on h_f. Such a term is conveniently represented by regarding the term $\frac{B}{F\left(q\right)}$ as a coefficient at each mesh point, and approximating $\frac{{\partial }^2{h}_f}{\partial x^2}$ by applying central differencing in space.

$F{\left(q\right)}_i^n=a-\frac{bK}{2\Delta x}\left(h_f{}_{i+1}^n-h_f{}_{i-1}^n\right)$ (29)

It should be noted that Equation (29) is expressed in dimensional form as the evaluation of the scalar function F(q) requires determining the flow rate through the fractures at all mesh points.

The derivative boundary conditions (25) and (26) are discretized by the backward finite difference formulation.

The discretized governing equations and boundary conditions yield a system of simultaneous algebraic equations for the blocks and the fractures, which is solved successively by backward substitution using the Thomas algorithm. Consequently, the dimensionless drawdowns in both the matrix blocks and the fractures are obtained as functions of dimensionless time and space. The entire numerical scheme is implemented in a FORTRAN code.

Table 1. Typical aquifer parameters.

Parameter
Width of aquifer, L	800 m
Thickness of aquifer, B	60 m
Initial piezometric head, h0	65 m
Constant drawdown in stream, s0	2 m
Transmissivity of fractures, Tf	280 m2/day
Transmissivity of blocks, Tb	1.75 × 10−2 m2/day
Storage coefficient of fractures, Sf	1.4 × 10−7
Storage coefficient of blocks, Sb	1.4 × 10−6
Fluid transfer rate, ε	6.25*10−3 m−2

The dimensionless drawdowns in both the fractures and the blocks are converted into dimensional drawdowns to enable the computation of flow rates at all mesh points to evaluate the scalar function F(q).

The calculations are conducted using a range of aquifer parameters which are given by [3]. Typical aquifer parameters are given in Table 1.

Furthermore, the following dimensionless parameters are introduced to characterize the aquifer and to facilitate a clearer understanding of the influence of the aquifer parameters on flow behavior.

$\eta =\frac{S_b}{S_f}$ (30)

where η is storativity contrast;

$\kappa =\frac{T_b}{T_f}$ (31)

where κ is conductivity contrast; and

$\delta =\frac{4\epsilon T_b{S}_f{L}^2}{T_f{S}_b}=4\epsilon L^2\frac{\kappa }{\eta }=4\epsilon L^2\frac{{\nu }_b}{{\nu }_f}$ (32)

where δ is diffusivity contrast; ν_b and ν_f are hydraulic diffusivities for the blocks and fractures respectively.

3.2. Von Neumann Stability Analysis

This section focuses exclusively on the von Neumann Stability Analysis, as the Discrete Perturbation Analysis has not provided reliable results and a conclusive explanation of the observed behavior.

Von Neumann stability analysis is a commonly used procedure for determining the stability requirements of finite difference equations. Practical problems typically involve variable coefficients, non-linearities and complicated types of boundary conditions. In these cases, the method can only be applied locally and with the non-linearities temporarily frozen. In the von Neumann method, the errors distributed along grid lines at one time level are expanded as Fourier series. Then the stability or instability of the computational algorithm is determined by considering whether separate Fourier components of the error distribution decay or amplify in progressing to the next time level [15].

Stability of the finite difference equations for the fractures are analyzed by von Neumann method. The Fourier component for $u_i^n$ is assumed as:

$u_i^n=U^n{\text{e}}^{IP\left(\Delta x\right)i}$ (33)

where $i=\sqrt{-1}$ , $U^n$ is the amplitude at time level n; P is the wave number in the x-direction; i.e. $k_x=\frac{2\pi }{P}$ , where $k_x$ is the wave length. Similarly:

$u_i^{n+1}=U^{n+1}{\text{e}}^{IP\left(\Delta x\right)i}$ (34)

and

$u_{i\pm 1}^n=U^n{\text{e}}^{IP\left(\Delta x\right)\left(i\pm 1\right)}$ (35)

If the phase angle is defined as

$\phi =P\Delta x$ (36)

Then,

$u_i^n=U^n{\text{e}}^{I\phi i}$ (37)

$u_i^{n+1}=U^{n+1}{\text{e}}^{I\phi i}$ (38)

and

$u_{i\pm 1}^n=U^n{\text{e}}^{I\phi \left(i\pm 1\right)}$ (39)

The solution of the finite difference equation is expanded in a Fourier series. After some mathematical operations, the finite difference equations have the form:

$U^{n+1}=G{U}^n$ (40)

where G is the amplification factor.

The decay or growth of the amplification factor indicates whether the numerical algorithm is stable. The necessary condition for stability is:

$\left|G\right|\le 1$ (41)

This condition is necessary and sufficient for two time level difference equations, but it is not always sufficient for three or more level equations, although it is always necessary.

Fourier components for the dimensionless piezometric heads for the fractures, i.e. $z_{f_i^n}=U^n{\text{e}}^{I\phi i}$ , $z_{f_i^{n+1}}=U^{n+1}{\text{e}}^{I\phi i}$ and $z_{f_{i\pm 1}^n}=U^n{\text{e}}^{I\phi \left(i\pm 1\right)}$ , are inserted in the finite difference equation for the fractures assuming that F(q) is constant.

Indeed F(q) = f(x, t) not constant and imposes non-linearity in the partial differential equation for the fractures. Reference [15] stated that von Neumann method gives some guidance for stability if the non-linearity is frozen locally. Furthermore, numerical calculations indicate that the value of F(q) does not change considerably. Therefore, F(q) can be assumed as a constant in order to be able to apply von Neumann method.

$\begin{array}{l}{U}^{n+1}{\text{e}}^{I\phi i}-U^n{\text{e}}^{I\phi i}\\ =\frac{B}{F\left(q\right)}\frac{1}{T_f}\frac{\Delta \theta }{\Delta {\lambda }^2}[\frac{1}{2}\left(U^{n+1}{\text{e}}^{I\phi i+1}-2U^{n+1}{\text{e}}^{I\phi i}+U^{n+1}{\text{e}}^{I\phi i-1}\right)\\ +\frac{1}{2}\left(U^n{\text{e}}^{I\phi i+1}-2U^n{\text{e}}^{I\phi i}+U^n{\text{e}}^{I\phi i-1}\right)]-\frac{T_b}{T_f}\Delta \theta L^2\epsilon U^n{\text{e}}^{I\phi i}\end{array}$ (42)

In order to obtain the amplification factor G, the following mathematical operations are carried out:

$\begin{array}{l}{U}^{n+1}-U^n\\ =\frac{B}{F\left(q\right)}\frac{1}{T_f}\frac{\Delta \theta }{\Delta {\lambda }^2}[\frac{1}{2}\left(U^{n+1}{\text{e}}^{I\phi }-2U^{n+1}+U^{n+1}{\text{e}}^{-I\phi }\right)\\ +\frac{1}{2}\left(U^n{\text{e}}^{I\phi }-2U^n+U^n{\text{e}}^{-I\phi }\right)]-\frac{T_b}{T_f}\Delta \theta L^2\epsilon U^n\end{array}$ (43)

Notice that

$\cos \phi =\frac{{\text{e}}^{I\phi }+{\text{e}}^{-I\phi }}{2}$ (44)

Then,

$\begin{array}{l}{U}^{n+1}-U^n\\ =\frac{B}{F\left(q\right)}\frac{1}{T_f}\frac{\Delta \theta }{\Delta {\lambda }^2}\left[\frac{1}{2}\left(2U^{n+1}\cos \phi -2U^{n+1}\right)+\frac{1}{2}\left(2U^n\cos \phi -2U^n\right)\right]\\ -\frac{T_b}{T_f}\Delta \theta L^2\epsilon U^n\end{array}$ (45)

$\begin{array}{l}{U}^{n+1}-U^n\\ =\frac{B}{F\left(q\right)}\frac{1}{T_f}\frac{\Delta \theta }{\Delta {\lambda }^2}\left[\left(U^{n+1}\cos \phi -U^{n+1}\right)+\left(U^n\cos \phi -U^n\right)\right]-\frac{T_b}{T_f}\Delta \theta L^2\epsilon U^n\end{array}$ (46)

$\begin{array}{l}{U}^{n+1}-U^{n+1}\frac{B}{F\left(q\right)}\frac{1}{T_f}\frac{\Delta \theta }{\Delta {\lambda }^2}\left(\cos \phi -1\right)\\ =U^n+U^n\frac{B}{F\left(q\right)}\frac{1}{T_f}\frac{\Delta \theta }{\Delta {\lambda }^2}\left(\cos \phi -1\right)-\frac{T_b}{T_f}\Delta \theta L^2\epsilon U^n\end{array}$ (47)

$\begin{array}{l}{U}^{n+1}\left[1-\frac{B}{F\left(q\right)}\frac{1}{T_f}\frac{\Delta \theta }{\Delta {\lambda }^2}\left(\cos \phi -1\right)\right]\\ =U^n\left[1+\frac{B}{F\left(q\right)}\frac{1}{T_f}\frac{\Delta \theta }{\Delta {\lambda }^2}\left(\cos \phi -1\right)-\frac{T_b}{T_f}\Delta \theta L^2\epsilon \right]\end{array}$ (48)

Recall that

$U^{n+1}=G{U}^n$ (49)

Therefore,

$G=\frac{\left[1+\frac{B}{F\left(q\right)}\frac{1}{T_f}\frac{\Delta \theta }{\Delta {\lambda }^2}\left(\cos \phi -1\right)-\frac{T_b}{T_f}\Delta \theta L^2\epsilon \right]}{\left[1-\frac{B}{F\left(q\right)}\frac{1}{T_f}\frac{\Delta \theta }{\Delta {\lambda }^2}\left(\cos \phi -1\right)\right]}$ (50)

$G=\frac{\left[1-\frac{B}{F\left(q\right)}\frac{1}{T_f}\frac{\Delta \theta }{\Delta {\lambda }^2}\left(1-\cos \phi \right)-\frac{T_b}{T_f}\Delta \theta L^2\epsilon \right]}{\left[1+\frac{B}{F\left(q\right)}\frac{1}{T_f}\frac{\Delta \theta }{\Delta {\lambda }^2}\left(1-\cos \phi \right)\right]}$ (51)

Since $\frac{T_b}{T_f}$ and $\Delta \theta$ are very small, the term $\frac{T_b}{T_f}\Delta \theta L^2\epsilon$ is always much less than 1. Consequently, the absolute value of Equation (51) remains below 1, ensuring that the finite difference equation for the fractures is unconditionally stable.

Fourier components for the dimensionless piezometric heads for the blocks, i.e. $z_{b_i^n}=U^n{\text{e}}^{I\phi i}$ , $z_{b_i^{n+1}}=U^{n+1}{\text{e}}^{I\phi i}$ and $z_{b_{i\pm 1}^n}=U^n{\text{e}}^{I\phi \left(i\pm 1\right)}$ , are inserted in the finite difference equation for the blocks, and the same procedure is applied to obtain the amplification factor, G.

$G=\frac{\left[1-\frac{S_f}{S_b}\frac{T_b}{T_f}\frac{\Delta \theta }{\Delta {\lambda }^2}\left(1-\cos \phi \right)-\frac{S_f}{S_b}\frac{T_b}{T_f}\Delta \theta L^2\epsilon \right]}{\left[1+\frac{S_f}{S_b}\frac{T_b}{T_f}\frac{\Delta \theta }{\Delta {\lambda }^2}\left(1-\cos \phi \right)\right]}$ (52)

4. Results and Discussion

Simulations are performed for various combinations of storativity, hydraulic conductivity and diffusivity contrasts, as well as for different combinations of time step and spatial step sizes, which yield different values of α, where

$\alpha =\frac{\Delta \theta }{\Delta {\lambda }^2}$ (53)

Firstly, $\Delta \theta$ is fixed to 1 × 10⁻⁵ and α is changed by changing $\Delta \lambda$ . Dimensionless buildup (drawdown) versus dimensionless time curves were generated for the fractures. These curves were compared with the curves available from the analytical solution presented by [3] at a specified dimensionless location $\lambda =0.2$ .

Figure 2 shows Case 1, for which η = 10, δ = 5 and κ = 3.125 × 10⁻⁵ at λ = 0.2. As evident from the figure that numerical and analytical solutions give similar results when α = 2.09. For this particular value of α, there is a slight difference between the curves of numerical and analytical solutions at early times, but the curves perfectly overlap at later times. The other values of α give rather different results and their curves do not fit the analytical data.

Figure 2. Z_f versus θ curves at λ = 0.2 for different α values (Case 1).

Figure 3 shows Case 2, for which η = 10, δ = 10 and κ = 6.25 × 10⁻⁵ at λ = 0.2. The curve of α = 2.09 again fits best the curve of the analytical data at λ = 0.2. However, for this case, the difference between these curves is more significant at early times, and the numerical solution gives a small delayed response. On the other hand, the numerical and analytical curves perfectly overlap at later times. The curves for the values of α other than 2.09 are separated from the analytical curve.

Figure 3. Z_f versus θ curves at λ = 0.2 for different α values (Case 2).

Figure 4. Z_f versus θ curves at λ = 0.2 for different α values (Case 3).

Figure 4 demonstrates Case 3, for which η = 10, δ = 1 and κ = 6.25 × 10⁻⁵ at λ = 0.2. The curve of α = 2.09 is the most similar curve to the curve of analytical solution. However, there is significant difference between these curves except for the very late times. On the other hand, they follow similar patterns.

Next, Δλ is fixed to 0.0022 and α is changed by altering Δθ. Figure 5 shows Case 4, for which η = 10, δ = 10 and κ = 6.25 × 10⁻⁵ at λ = 0.2. It can be concluded from the figure that the numerical and analytical solutions yield similar results when α = 2.09. The difference between the curve of α = 2.09 and the curve of the analytical data is more significant at early times, and the numerical solution gives a small delayed response.

Figure 5. Z_f versus θ curves at λ = 0.2 for different α values (Case 4).

Figure 6. Z_f versus θ curves at λ = 0.4 for different α values (Case 5).

Figure 6 illustrates Case 5, for which η = 10, δ = 10 and κ = 6.25 × 10⁻⁵, corresponding to a different dimensionless location, λ = 0.4. Once again, the numerical and analytical solutions exhibit the closest agreement when α = 2.09. Furthermore, the deviation between the numerical and analytical results at early times is smaller at this dimensionless location because the hydraulic gradient is relatively low and, consequently, the nonlinear effect becomes less significant.

Figure 7 illustrates Case 6, for which η = 10, δ = 5 and κ = 3.125 × 10⁻⁵, corresponding to a different dimensionless location, λ = 0.8. The numerical and analytical results are in very close agreement, with an almost perfect match observed at α = 2.09. Moreover, the difference between the numerical and analytical curves is very small even at early times.

The comparison indicates that the numerical solution exhibits the best visual agreement with the analytical results when α = 2.09. However, it would be necessary to prove it quantitatively. To this end, root mean square (RMS) error for six values α were calculated for five cases employing Equation (54) [16].

Figure 7. Z_f versus θ curves at λ = 0.8 for different α values (Case 6).

$\text{RMS error}={\left[\frac{1}{N}{\displaystyle \sum }_1^N{\left|z_{fanalytical}-z_{fnumerical}\right|}^2\right]}^{\frac{1}{2}}$ (54)

where N is the number of points compared.

Table 2 presents the RMS errors computed at the points for which analytical solutions are available. The results show that, for all five cases considered, the minimum RMS error is achieved when α = 2.09. This indicates that the discrepancy between the numerical and analytical solutions is minimized at this value of α. These findings further confirm the excellent agreement between the numerical and analytical results obtained for Case 6.

The amplification factor, G, associated with the finite difference equations for the fractures, has been evaluated for the six cases presented above over a range of α values. The resulting G versus α curves were plotted to investigate the relationship between these parameters.

It is evident from Figure 8 that the amplification factor curves for all cases attain a minimum value when α is approximately 2, corresponding to the closest agreement between the numerical and analytical solutions, regardless of whether Δθ or Δλ is held constant while the other is varied.

For Cases 1, 2 and 3, where λ = 0.2 and Δθ is held constant at 0.00001, G versus α curves show only slight variations, with their minima close to zero.

For Cases 5 and 6, where λ = 0.4 and λ = 0.8 respectively and Δθ is held constant at 0.00001, G versus α curves exhibit noticeable divergence for α > 10. In both cases, the curves attain their minimum values at approximately G = 0.07.

For Case 4, where Δλ is fixed at 0.0022, the curve differs to some extend and its minimum occurs around G = 0.1.

These findings quantitatively validate the selection of α = 2.09 as the value that produces the best agreement between the numerical and analytical solutions.

It is worth noting that a minimum amplification factor may indicate that the numerical scheme neither introduces excessive damping nor causes artificial amplification of the underlying physical modes, particularly in response to the initial discontinuity that is the step rise in piezometric head imposed at the stream boundary. Consequently, the discrete solution preserves the amplitudes of the dominant Fourier components, and provides a more accurate representation of the true physical behavior.

The same procedure is applied to the finite difference equations for the matrix blocks. For all cases, the amplification factor, G, is found to be very close to 1. Therefore, no clear explanation is obtained for why the numerical solution achieves optimal accuracy at α = 2.09.

Table 2. RMS errors.

RMS error	α = 0.028	α = 0.5	α = 1	α = 2.09	α = 6.4	α = 25.6
Case 1	0.178	0.112	0.064	0.012	0.145	0.332
Case 2	0.314	0.236	0.164	0.102	0.136	0.195
Case 3	0.170	0.141	0.120	0.088	0.157	0.212
Case 5	0.298	0.185	0.162	0.058	0.148	0.176
Case 6	0.308	0.243	0.192	0.009	0.152	-

Figure 8. G versus α curves for all cases.

This behavior may be attributed to the initially slow buildup response of the matrix blocks, which act as large storage reservoirs that gradually absorb and release water. As a result, pressure equilibration occurs more slowly and the transient modes decay at a reduced rate. Accordingly, the matrix blocks respond more slowly to the increase in piezometric head, limiting the ability of the discrete solution to accurately preserve the amplitudes of the dominant Fourier modes.

In addition to this, the Crank-Nicolson implicit scheme is unconditionally stable, but not strongly dissipative. For certain modes, especially when α is large, $\left|G\right|\approx 1$ , which means that the oscillatory components are propagated as weakly damped oscillations in time [17]. This behavior is likely reflected in the finite difference equations governing the matrix blocks. Since the computed amplification factor is very close to unity, numerical perturbations experience only weak damping and may therefore persist within the matrix domain for relatively long periods.

5. Conclusions

The validity of the proposed numerical scheme is assessed. The results obtained are compared in detail with an established analytical solution, thereby verifying accuracy. In parallel, von Neumann stability analysis is conducted within a systematic quantitative framework to investigate the stability behavior and to elucidate the scheme’s response to different governing parameters.

The comparison of the analytical and numerical results for six different cases indicates that the analytical and numerical solutions agree closely when α = 2.09.

To further strengthen the validation of the numerical model, RMS errors are computed at the locations for which analytical solutions are available. For five cases examined, the smallest RMS error is consistently obtained when α = 2.09, indicating that this value yields the closest agreement between the analytical and numerical solutions.

The magnitude of the amplification factor, G, associated with the finite difference equations for the fractures has been computed for the six cases considered over a range of α values. The resulting G versus α curves are examined to assess the influence of α on the amplification behavior of the scheme. For all cases, the amplification factor reaches its minimum at α = 2.09, which also corresponds to the minimum discrepancy between the analytical and numerical solutions.

The combined strategy ensures a thorough evaluation and confirmation of the numerical solution’s reliability and accuracy. It also offers a comprehensive basis for assessing the performance of the scheme under varying governing parameters.

Conflicts of Interest

The author declares no conflicts of interest regarding the publication of this paper.

References

[1]	Altinörs, A. and Önder, H. (2008) A Double-Porosity Model for a Fractured Aquifer with Non-Darcian Flow in Fractures. Hydrological Sciences Journal, 53, 868-882.[CrossRef]
[2]	Altınörs, A.A. (2005) Non-Darcian Flow in a Fractured Aquifer. Ph.D. Thesis, Middle East Technical University.
[3]	Önder, H. (1998) One-Dimensional Transient Flow in a Finite Fractured Aquifer System. Hydrological Sciences Journal, 43, 243-265.[CrossRef]
[4]	Bear, J. (1979) Hydraulics of Groundwater. McGraw-Hill.
[5]	Bai, M., Elsworth, D. and Roegiers, J. (1993) Multiporosity/Multipermeability Approach to the Simulation of Naturally Fractured Reservoirs. Water Resources Research, 29, 1621-1633.[CrossRef]
[6]	Choi, E.S., Cheema, T. and Islam, M.R. (1997) A New Dual-Porosity/Dual-Permeability Model with Non-Darcian Flow through Fractures. Journal of Petroleum Science and Engineering, 17, 331-344.[CrossRef]
[7]	Huyakorn, P.S. and Pinder, G.F. (1983) Computational Methods in Subsurface Flow. Academic Press.
[8]	Wu, Y.S. (2002) Numerical Simulation of Single-Phase and Multiphase Non-Darcy Flow in Porous and Fractured Reservoirs. Transport in Porous Media, 49, 209-240.[CrossRef]
[9]	Barenblatt, G.I., Zheltov, I.P. and Kochina, I.N. (1960) Basic Concepts in the Theory of Seepage of Homogeneous Liquids in Fissured Rocks [Strata]. Journal of Applied Mathematics and Mechanics, 24, 1286-1303.[CrossRef]
[10]	Streltsova, T.D. (1975) Unsteady Unconfined Flow into a Surface Reservoir. Journal of Hydrology, 27, 95-110.[CrossRef]
[11]	Warren, J.E. and Root, P.J. (1963) The Behavior of Naturally Fractured Reservoirs. Society of Petroleum Engineers Journal, 3, 245-255.[CrossRef]
[12]	Hirsch, C. (1989) Numerical Computation of Internal and External Flows. Vol. 1, John Wiley and Sons Interscience Publications.
[13]	Fadugba, E.F., Edognabya, H.O. and Zelibe, S.C. (2013) Crank-Nicolson Method for Solving Parabolic Partial Differential Equations. International Journal of Applied Mathematics and Modeling, 1, 8-23.
[14]	Morton, K.W. and Mayers, D.F. (2005) Numerical Solution of Partial Differential Equations. 2nd Edition, Cambridge University Press.[CrossRef]
[15]	Fletcher, C.A.J. (1991) Computational Techniques for Fluid Dynamics I. Springer-Verlag.
[16]	Mathews, J.H. (1988) Numerical Methods for Computer Science, Engineering and Mathematics. Prentice-Hall International Editions.
[17]	Britz, D., Østerby, O. and Strutwolf, J. (2003) Damping of Crank-Nicolson Error Oscillations. Computational Biology and Chemistry, 27, 253-263.[CrossRef] [PubMed]