Keywords

1 Introduction

The flow and solute transport in fractured porous media is very important in many applications such as contaminants migration in fractured aquifer systems. The fractured porous medium consists of two domains, namely, matrix blocks and fractures. The fractures are more permeable than the matrix blocks, but they contain very little fluid than the matrix. Different scales are therefore invoked in fractured porous media such as discrete fracture models (DFMs) and dual continuum models, can describe flow and transport in fractured porous media. The model of solute transport in fractured and porous media has been solved analytical by Park and Lee [1] and Roubinet et al. [2]. In order to model solute transport in fractured porous media, Bodin et al. [3] and Graf and Therrien [4] have used the discrete network model; and Refs. [5,6,7] have used equivalent continuum models, while Refs. [8, 9] have used the continuum model. In order to represent the fractures explicitly in the fractured porous media the discrete-fracture model (DFM) has been used. This procedures can remove the contrast of the length-scale resulting of the direct representation of the fracture aperture as in the dual-porosity model.

Time discretization, on the other hand, has a significant impact on the efficiency of numerical solutions. Therefore, the use of traditional single-scale time schemes is limited by the rapid variation in pressure and concentration in matrix or fracture, where applicable. The multi-scale time-splitting technique is therefore considered to be one of the major improvements in the treatment of the gap between pressure and concentration. In a number of publications, such as [10,11,12,13], the multi-scale time splitting method was considered. El-Amin et al. [14] have developed a discrete-fracture-model with multi-scale time-splitting two-phase flow including nanoparticles transport in fractured porous media. In this work, we develop a discrete-fracture-model with multi-scale time-splitting of solute transport in fractured porous media. The modeling and mathematical formulation is considered in the second section. In Sect. 3, the time-stepping technique with spatial discretization have been presented. The fourth section is devoted to numerical test, and then, the conclusions are given in the last section.

2 Modeling and Formulation

2.1 Governing Equations

This paper considers the problem of mass transfer and flow in fractured porous media. The system of equations consists of continuity, momentum, and concentration.

Momentum Conservation (Darcy’s Law):

$$\begin{aligned} \mathbf {u} = -\mathbf {K} \nabla \varPhi \end{aligned}$$
(1)

where \(\mathbf {K}\) is the permeability tensor \(\mathbf {K}=\frac{k}{\mu }\mathbf {I}\), \(\mathbf {I}\) is the identity matrix and \(k/\mu \) is a positive real number. \(\mathbf {u}, \varPhi , \mu , k\) are, respectively, the velocity, the pressure, the viscosity and the permeability.

Mass Conservation:

$$\begin{aligned} \nabla \cdot \mathbf {u} = q, \end{aligned}$$
(2)

where q is the external mass flow rate.

Mass Transfer: The solute transport equation in porous media may be given as,

$$\begin{aligned} \phi \frac{\partial C}{\partial t} + \nabla \cdot \left( \mathbf {u} C - \phi D \nabla C \right) = Q_c, \end{aligned}$$
(3)

where C is the solute concentrations, \(\phi \) is the porosity, D is the diffusion coefficient, and \(Q_c\) is the rate of change of volume belonging to a source/sink term.

2.2 Initial and Boundary Conditions

Consider the computational domain \(\varOmega \) with the boundary \(\partial \varOmega \) which is subjected to Dirichlet \(\varGamma _D\) and Neumann \(\varGamma _N\) boundaries, where \(\partial \varOmega = \varGamma _D \cup \varGamma _N\) and \(\varGamma _D \cap \varGamma _N=\) ø. At the beginning of the injection process, we have,

$$\begin{aligned} C = 0 \quad \mathrm{in} \quad \varOmega \quad \quad \mathrm{at} \quad t=0, \end{aligned}$$
(4)

The boundary conditions are given as,

$$\begin{aligned} P = P^D \quad \mathrm{on} \quad \varGamma _D, \end{aligned}$$
(5)
$$\begin{aligned} \mathbf {u} \cdot \mathbf {n} = q^N, \quad C = C_0, \quad \mathrm{on} \quad \varGamma _N. \end{aligned}$$
(6)

where \(\mathbf {n}\) is the outward unit normal vector to \(\partial \varOmega \), \(P^D\) is the pressure on \(\varGamma _D\) and \(q^N\) the imposed inflow rate on \(\varGamma _N\), respectively.

2.3 Discrete Fracture Model

In the discrete-fracture-model (DFM), the fracture gridcells are simplified to represent as the interfaces of the matrix gridcell and fractures are surrounded by matrix blocks. Thus, the dimension of fracture reduced by one than the dimension of matrix, i.e., if matrix is of n-dimension, then, fracture is of (\(n-1\))-dimension. The domain is decomposed into the matrix domain, \(\varOmega _{\text {m}}\) and fracture domain, \(\varOmega _{\text {f}}\). The pressure equation in the matrix domain is given by,

$$\begin{aligned} -\nabla \cdot \mathbf{K} _{\text {m}}\nabla \varPhi _{\text {m}} = q_m, \end{aligned}$$
(7)

Assuming that the pressure along the fracture width are constants, and by integration, the pressure equation in the fracture becomes,

$$\begin{aligned} -\nabla \cdot \mathbf{K} _{\text {f}}\nabla \varPhi _{\text {f}} = q_{\text {f}}+Q_{\text {f}}, \end{aligned}$$
(8)

The matrix-fracture interface condition is given by,

$$\begin{aligned} \varPhi _{\text {m}}= \varPhi _{\text {f}}, ~~~\text {on} ~~\partial \varOmega _{\text {m}}\cap \varOmega _{\text {f}}. \end{aligned}$$
(9)

where the subscript m represents the matrix domain, while the subscript f represents the fracture domain. \(Q_{\text {f}}\) is the mass transfer across the matrix-fracture interfaces.

The solute transport equation in the matrix domain may be expressed as,

$$\begin{aligned} \phi _\text {m} \frac{\partial C_\text {m} }{\partial t} + \nabla \cdot \left( \mathbf {u}_{\text {m}} C_\text {m} - \phi _\text {m} D\nabla C_\text {m} \right) = Q_{c,\text {m}}, \end{aligned}$$
(10)

where \(C_\text {m}\) is the concentration in the matrix domain. \(Q_{c,\text {m}}\) is the rate of change of volume belonging to a source/sink term in the matrix domain. On the other hand, the solute transport equation in fractures is represented by,

$$\begin{aligned} \phi _\text {f} \frac{\partial C_\text {f} }{\partial t} + \nabla \cdot \left( \mathbf {u}_{\text {f}} C_\text {f} - \phi _\text {f} D\nabla C_\text {f} \right) = Q_{c,\text {f}} + Q_{c,\text {f}}, \end{aligned}$$
(11)

where \(C_\text {f}\) is the mass concentration in the fracture domain. \(Q_{c,\text {f}}\) is the rate of change of particle belonging to a source/sink term in the fracture domain. \(Q_{c,\text {f}}\) represents the rate of change of volume across the matrix-fracture interfaces. The interface condition of the solute concentration is,

$$\begin{aligned} C_{\text {m}}= C_{\text {f}}, ~~~\text {on} ~~~\partial \varOmega _{\text {m}}\cap \varOmega _{\text {f}}. \end{aligned}$$
(12)

3 Multiscale Time-Splitting and Spatial Discretization

In the multiscale time-splitting method, we employ a different time step-size for each time derivative term as they have different physics. For example, the time-step size for the pressure can be larger than it of the solute concentration. Also, the fractures pressure may be has a larger time-step size than one for the matrix. So, we may use a small time step-size for the pressure in fractures, and so on. We use the CCFD method for the spatial discretization. The CCFD method is locally conservative and equivalent to the quadratic mixed finite element method.

3.1 Multiscale Time-Splitting Approach for Pressure

Now, let us introduce the time discretization for the pressure in the matrix domain. The total time interval [0, T] is divided into \(N_{p,\text {m}}\) steps, i.e., \(0=t^0<t^1<\cdots <t^{N_{p,\text {m}}}=T\) and the time step length is \(\triangle t^i=t^{i+1}-t^i\). Therefore, we divide each subinterval \((t^i, t^{i+1}]\) into \(N_{p,\text {f}}\) sub-subintervals as \((t^i,t^{i+1}]=\bigcup _{j=0}^{N_{p,\text {f}}-1}(t^{i,j},t^{i,j+1}]\), where \(t^{i,0}=t^{i}\) and \(t^{i,N_{p,\text {f}}}=t^{i+1}\) and \(\triangle t^{i,j}=t^{i,j+1}-t^{i,j}\). In the following, b refers to the boundary of the matrix gridcells K such that its area is |K|, and |b| is its length. \(\mathbf{n} _{b}\) is a unit normal vector pointing from K to \(K'\) on each interface \(b\in \partial K\cap \partial K'\). The flux across the boundary b of the gridcell K is denoted by \(\mathbf {\xi }\). \(d_{K,b}\) is the distance from the central points of the cell K and the cell boundary b. \(d_{K,K'}\) is distance between the central points of the cells K and \(K'\). When b is located on the entire domain boundary, the pressure is provided by Dirichlet boundary conditions, \(b\in \varGamma ^D\). Otherwise, the Neumann conditions \(b\in \partial \varOmega ^{N}\) is used to calculate fluxes.

The pressure equation in the matrix domain and fractures is written, respectively as,

$$\begin{aligned} -\nabla \cdot \mathbf{K} _{\text {m}}\nabla \varPhi _{\text {m}}^{i+1} = q_{\text {m}}^{i+1}, \end{aligned}$$
(13)

and

$$\begin{aligned} -\nabla \cdot \mathbf{K} _{\text {f}}\nabla \varPhi _{\text {f}}^{i+1} = q_{\text {f}}^{i+1}+Q_{\text {f}}^{i+1}. \end{aligned}$$
(14)

Now, applying the CCFD scheme on (13), one obtains,

$$\begin{aligned} \sum _{b\in \partial K} \mathbf{{\xi }_{a,\mathrm{{m}},b}^{i+1} = q_{\mathrm{{m}},K}^{i+1}|K|}, \end{aligned}$$
(15)

If \(b\in \partial K\cap \partial K'\) and \(b\nsubseteq \varOmega _{\text {f}}\), the fluxes in (15) are given by,

$$\begin{aligned} \mathbf{{\xi }_{a,\mathrm{{m}},b}^{i+1}=-|b|\chi _{t,b}^i \frac{ \Phi _{\mathrm{{m}},K'}^{i+1}- \Phi _{\mathrm{{m}},K}^{i+1}}{d_{K,K'}}}, \end{aligned}$$
(16)

where \(\chi _{t,b}^i\) is given by

$$\begin{aligned} \chi _{t,b}^i=\frac{d_{K,K'}{} \mathbf{K} _{\text {m},K}{} \mathbf{K} _{\text {m},K'}}{d_{K,b}{} \mathbf{K} _{\text {m},K'}+d_{K',b}{} \mathbf{K} _{\text {m},K}}, \end{aligned}$$
(17)

On the other hand, if \(b\in \varOmega _{\text {f}}\cap \partial K\cap \partial K'\) and b is a gridcell of the fracture system, we have,

$$\begin{aligned} \mathbf{{\xi }_{a,\mathrm{{m}},b}^{i+1}\equiv {\xi }_{a,\mathrm{{m}},b,K}^{i+1}=-|b|\chi _{t,\mathrm{{mf}},b}^i \frac{ \Phi _{\mathrm{{f}},b}^{i+1}- \Phi _{\mathrm{{m}},K}^{i+1}}{d_{K,b}+\frac{\epsilon }{2}}}, \end{aligned}$$
(18)

where \(\chi _{t,\text {mf},b}^i\) is given by,

$$\begin{aligned} \chi _{t,\text {mf},b}^i=\frac{(d_{K,b}+\frac{\epsilon }{2})\mathbf{K} _{\text {f},K} \mathbf{K} _{\text {m},K}}{\frac{\epsilon }{2} \mathbf{K} _{\text {m},K}+d_{K,b} \mathbf{K} _{\text {f},K}}, \end{aligned}$$
(19)

Similarly, let b be a gridcell of the fracture network. Equation (14) may be discretized by the CCFD method to get,

$$\begin{aligned} \sum _{\gamma \in \partial b} \mathbf{{\xi }_{a,\mathrm{{f}},\gamma }^{i+1} = q_{\mathrm{{f}},b}^{i+1}|b| + Q_{\mathrm{{f}},b}^{i+1}|b|}, \end{aligned}$$
(20)

where \(\gamma \) is the face of the gridcell b in the fracture network. \(\mathbf {\xi }\) is the flux across the boundary \(\gamma \) of the fracture gridcell b. The matrix-fracture transfer is treated as a source term in the fracture system.

$$\begin{aligned} Q_{\text {f},b}^{i+1}= -(Q_{\text {f},b,K}^{i+1}+Q_{\text {f},b,K'}^{i+1})/\epsilon , \end{aligned}$$
(21)
$$\begin{aligned} Q_{\text {f},b,K}^{i+1}=\mathbf{{\xi }_{a,\mathrm{{m}},b,K}^{i+1}}, \end{aligned}$$
(22)
$$\begin{aligned} Q_{\text {f},b,K'}^{i+1}=\mathbf{{\xi }_{a,\mathrm{{m}},b,K'}^{i+1}}, \end{aligned}$$
(23)

where \(\mathbf{{\xi }_{a,\mathrm{{m}},b}^{i+1}}\) is defined in (18).

In the case of multiple fractures that connected by the interface \(\gamma \). Assume that \(\varLambda _\gamma \) is the set of the fracture grid cells joint by \(\gamma \). The mass conservation equation discretization is,

$$\begin{aligned} \sum _{b\in \varLambda _\gamma } \mathbf{{\xi }_{a,\mathrm{{f}},\gamma ,e}^{i+1} = 0}, \end{aligned}$$
(24)

where \(\mathbf{{\xi }_{a,\mathrm{{f}},\gamma ,b}={\xi }_{a,\mathrm{{f}},\gamma }|_{\gamma \in b}}\) and \(\mathbf{{\xi }_{c,\mathrm{{f}},\gamma ,b}={\xi }_{c,\mathrm{{f}},\gamma }|_{\gamma \in b}}\).

The discretization of the total mass conservation of the matrix domain and the fractures network is represented as,

$$\begin{aligned}&\left[ \begin{array}{cc} \mathbf{A} _{a,\text {m,m}}^i &{} \mathbf{A} _{a,\text {m,f}}^{i} \\ \mathbf{A} _{a,\text {f,m}}^{i} &{} \mathbf{A} _{a,\text {f,f}}^{i} \\ \end{array} \right] \left[ \begin{array}{c} \varPhi _{\text {m}}^{i+1} \\ \varPhi _{\text {f}}^{i+1} \\ \end{array} \right] =\left[ \begin{array}{c} \mathbf{Q} _{ac,\text {m}}^{i+1} \\ \mathbf{Q} _{ac,\text {f}}^{i+1} \\ \end{array} \right] . \end{aligned}$$
(25)

The pressure equation in the fractures at each time subtime-step is given by,

$$\begin{aligned} -\nabla \cdot \mathbf{K} _{\text {f}}\nabla \varPhi _{\text {f}}^{i,j+1} = q_{\text {f}}^{i,j+1}+Q_{\text {f}}^{i,j+1}. \end{aligned}$$
(26)

It is obtained by using the CCFD scheme to (26) that

$$\begin{aligned} \sum _{\gamma \in \partial e} \mathbf{{\xi }_{a,\mathrm{{f}},\gamma }^{i,j+1} = q_{\mathrm{{f}},e}^{i,j+1}|e| + Q_{\mathrm{{f}},e}^{i,j+1}|e|}, \end{aligned}$$
(27)

The matrix-fracture transfer is given by,

$$\begin{aligned} Q_{\text {f},b}^{i,j+1} = -(Q_{\text {f},b,K}^{i,j+1} + Q_{\text {f},b,K'}^{i,j+1})/\epsilon , \end{aligned}$$
(28)
$$\begin{aligned} Q_{\text {f},b,K}^{i,j+1}=\mathbf{{\xi }_{a,\mathrm{{m}},e,K}^{i,j+1}}, \end{aligned}$$
(29)
$$\begin{aligned} Q_{\text {f},e,K'}^{i,j+1}=\mathbf{{\xi }_{a,\mathrm{{m}},b,K'}^{i,j+1}}, \end{aligned}$$
(30)

where

$$\begin{aligned} \mathbf{{\xi }_{a,\mathrm{{m}},b,K}^{i,j+1}\equiv {\xi }_{a,\mathrm{{m}},b}^{i,j+1}=-|e|\chi _{t,\mathrm{{mf}},b}^{i,j} \frac{ \Phi _{\mathrm{{f}},b}^{i,j+1}- \Phi _{\mathrm{{m}},K}^{i+1}}{d_{K,b}+\frac{\epsilon }{2}}}, \end{aligned}$$
(31)

along with,

$$\begin{aligned} \chi _{t,\text {mf},b}^{i,j}=\frac{(d_{K,b}+\frac{\epsilon }{2})\mathbf{K} _{\text {f},K}{} \mathbf{K} _{\text {m},K}}{\frac{\epsilon }{2}{} \mathbf{K} _{\text {m},K}+d_{K,b}{} \mathbf{K} _{\text {f},K}}, \end{aligned}$$
(32)

For the case of multiple fractures, the pressure equation may be given as,

$$\begin{aligned} \mathbf{A} _{\text {f}}^{i,j} \varPhi _{\text {f}}^{i,j+1} =\mathbf{Q} _{\text {f}}^{i,j}, \end{aligned}$$
(33)

where

$$\begin{aligned} \mathbf{A} _{\text {f}}^{i,j}=\mathbf{A} _{a,\text {f,f}}^{i,j}, \end{aligned}$$
(34)

and

$$\begin{aligned} \mathbf{Q} _{\text {f}}^{i,j}= & {} \mathbf{Q} _{ac,\text {f}}^{i,j+1} -\mathbf{A} _{a,\text {f,m}}^{i,j} \varPhi _{\text {m}}^{i+1}. \end{aligned}$$
(35)

At the time step \((t^i,t^{i+1})\), (25) is solved implicitly to get \( \varPhi _{\text {m}}^{i+1}\). Then, at the time substep \((t^{i,j},t^{i,j+1})\), we compute \( \varPhi _{\text {f}}^{i,j+1}\) using (33). After that, we calculate fluxes as explained below. For the boundary b of the matrix gridcell K, \(\mathbf{n} _{b}\) is the unit normal vector pointing towards outside K. If \(e\in \partial K\cap \partial K'\) and \(b\nsubseteq \varOmega _{\text {f}}\),

$$\begin{aligned} \mathbf{{\xi }_{a,\mathrm{{m}},b}^{i,j+1}=-|b|\chi _{t,b}^{i,j} \frac{ \Phi _{\mathrm{{m}},K'}^{i+1}- \Phi _{\mathrm{{m}},K}^{i+1}}{d_{K,K'}}}, \end{aligned}$$
(36)

where

$$\begin{aligned} \chi _{t,b}^{i,j}=\frac{d_{K,K'}{} \mathbf{K} _{\text {m},K}{} \mathbf{K} _{\text {m},K'}}{d_{K,b} \mathbf{K} _{\text {m},K'}+d_{K',b} \mathbf{K} _{\text {m},K}}. \end{aligned}$$
(37)

If \(b\in \varOmega _{\text {f}}\cap \partial K\cap \partial K'\) and b is a fracture gridcell.

3.2 Multiscale Time-Stepping of the Concentration Equation

On the other hand, as the solute concentration vary more rapidly than the pressures. We also use a smaller time-step size for the concentration in matrix domain and the smallest time-step size for the concentration in fractures. The backward Euler time discretization is used for the equation of concentration. Therefore, the system of governing equations is solved based on the multiscale time-splitting technique. Now, let us divide the time-step \((t^{i,j}, t^{i,j+1}]\) of the fractures pressure into \(N_{c,\text {m}}\) sub-steps such that \((t^{i,j}, t^{i,j+1}]=\bigcup _{k=0}^{N_{s,\text {m}}-1}(t^{i,j,k},t^{i,j,k+1}]\), \(t^{i,j,0}=t^{i,j}\) and \(t^{i,j,N_{s,\text {m}}}=t^{i,j+1}\). This time discretization is employed for the concentration in the matrix domain. Moreover, we use a smaller time-step size for the fracture concentration. Thus, we partition the time-step, \((t^{i,j,k}, t^{i,j,k+1}]\) into \(N_{c,\text {f}}\) time sub-steps as \((t^{i,j,k}, t^{i,j,k+1}]=\bigcup _{l=0}^{N_{c,\text {f}}-1}(t^{i,j,k,l},t^{i,j,k,l+1}]\), where \(t^{i,j,k,0}=t^{i,j,k}\) and \(t^{i,j,k,N_{c,\text {f}}}=t^{i,j,k+1}\). The concentration is computed implicitly as follow,

$$\begin{aligned} \phi _{\text {m}} \frac{C_\text {m}^{i,j,k+1} - C_\text {m}^{i,j,k}}{\varDelta t^{i,j,k}} + \nabla \cdot \left( \mathbf {u}_{\text {m}}^{i+1} C_\text {m}^{i,j,k+1} - \phi _{\text {m}} D \nabla C_\text {m}^{i,j,k+1} \right) = Q_{c,\text {m}}^{i,j,k+1} \end{aligned}$$
(38)

In a similar manner, we consider variation of the concentration in the fractures are faster than those in the matrix domain. So, the concentration in the fractures is expressed as follow,

$$\begin{aligned} \phi _{\text {f}} \frac{ C_\text {f}^{i,j,k,l+1} - C_\text {f}^{i,j,k,l}}{\varDelta t^{i,j,k,l}} + \nabla \cdot \left\{ \mathbf {u}_{\text {f}}^{i+1} C_\text {f}^{i,j,k,l+1} - \phi _{\text {f}} D \nabla C_\text {f}^{i,j,k,l+1} \right\} = Q_{c,\text {f}}^{i,j,k,l+1} + Q_{c,\text {m},\text {f}}^{i,j,k,l+1} \end{aligned}$$
(39)

We use the upwind CCFD method to discretize the concentration Eq. (38),

$$\begin{aligned} |K|\phi _{\text {m},K}\frac{ C_{\text {m},K}^{i,j,k+1} - C_{\text {m},K}^{i,j,k}}{\varDelta t^{i,j,k}} + \sum _{b\in \partial K} {\hat{C}}_{m,K}^{i,j,k+1} \mathbf {F}_{a,\text {m},b}^{i,j+1} + \sum _{b\in \partial K} \mathbf {F}_{D,\text {m},b}^{i,j,k+1} = Q_{c,\text {m},K}^{i,j,k+1}|K|. \end{aligned}$$
(40)

where

$$ \mathbf {F}_{\text {m},b}^{i,j+1} = \mathbf {u}_{\text {m},b}^{i,j+1} |b| $$

Let b be the interface between the matrix gridcells K and \(K'\); that is, \(b=\partial K\cap \partial K'\). If \(b\nsubseteq \varOmega _{\text {f}}\), the term \({\hat{C}}_{m,K}^{i,j,k+1}\) in (40) is given by,

$$\begin{aligned} \hat{C}_{m,K}^{i,j,k+1} =\left\{ \begin{array}{c} C_{\mathrm{{m}},K}^{i,j,k}, ~~~ \mathbf{{F}_{\mathrm{{m}},b}^{i,j+1}>0} , \\ C_{\mathrm{{m}},K'}^{i,j,k}, ~~~ \mathbf{{F}_{\mathrm{{m}},b}^{i,j+1}<0} .\\ \end{array} \right. \end{aligned}$$
(41)

Now for the diffusion term; if \(b\in \partial K\cap \partial K'\) and \(b\nsubseteq \varOmega _{\text {f}}\), the fluxes in (40) are given by,

$$\begin{aligned} \mathbf{F} _{D,\text {m},b}^{i,j,k+1} = -|b|\chi _{t,b}^{i,j,k+1} \frac{ C_{\text {m},K'}^{i,j,k+1}- C_{\text {m},K}^{i,j,k+1}}{d_{K,K'}}, \end{aligned}$$
(42)

where \(\chi _{t,b}^{i,j,k}\) is given by the harmonic mean as,

$$\begin{aligned} \chi _{t,b}^{i,j,k+1} = \frac{d_{K,K'} D_{\text {m},K}^{i,j,k} D_{\text {m},K'}^{i,j,k} \phi ^{i,j,k}_{\text {m},K} \phi ^{i,j,k}_{\text {m},K'} }{d_{K,b} D_{\text {m},K'}^{i,j,k} \phi ^{i,j,k}_{\text {m},K'} + d_{K',b}D_{\text {m},K}^{i,j,k} \phi ^{i,j,k}_{\text {m},K} }, \end{aligned}$$
(43)

On the other hand, if \(b\in \varOmega _{\text {f}}\cap \partial K\cap \partial K'\) and b is a gridcell of the fracture system, we have,

$$\begin{aligned} \mathbf{F} _{D,\text {m},b}^{i,j,k,l+1} \equiv \mathbf{F} _{D,\text {m},b,K}^{i,j,k+1} =-|b|\chi _{t,\text {mf},b}^{i,j,k} \frac{ C_{\text {f},b}^{i,j,k+1} - C_{\text {m},K}^{i,j,k+1}}{d_{K,b}+\frac{\epsilon }{2}}, \end{aligned}$$
(44)

where \(\chi _{t,\text {mf},b}^{i,j,k,l+1}\) is defined as,

$$\begin{aligned} \chi _{t,\text {mf},b}^{i,j,k,l+1} = \frac{(d_{K,b}+\frac{\epsilon }{2}) D_{\text {m},K}^{i,j,k} D_{\text {f},K}^{i,j,k} \phi ^{i,j,k}_{\text {m},K} \phi ^{i,j,k}_{\text {f},K} }{\frac{\epsilon }{2} D_{\text {m},K}^{i,j,k} \phi ^{i,j,k}_{\text {m},K} + d_{K,b} D_{\text {f},K}^{i,j,k} \phi ^{i,j,k}_{\text {f},K} }, \end{aligned}$$
(45)
Table 1. Physical and computational parameters

4 Numerical Tests

In order to examine the proposed scheme, three examples of fractured media with different dimensions, namely, 20 m \(\times \) 15 m \(\times \) 1 m, 10 m \(\times \) 10 m \(\times \) 1 m, and 10 m \(\times \) 10 m \(\times \) 1 m are presented with different multiple interconnected fractures as shown in Fig. 1. The matrix permeability is taken as 1 md in Examples (1) and (2), while in Example (3) is taken as 50 md. On the other hand, the fractures permeability is taken as \(10^6\) md for Examples (2) and (3), while in Example (1) is taken as \(10^5\) md. The total number of gridcells is 2500 in Examples (1) and (2), while in Example (3) is 3300. The solute is injected into a water aquifer with a rate of 0.1 PV with initial concentration of 0.1. The remaining of the physical and computational parameters are given for the three examples in Table 1. The outer (pressure) time-step size is taken as \(\varDelta t=1.9\), while we chose, \(N_{p,\text {f}}=5\), \(N_{c,\text {m}}=2\) and \(N_{c,\text {f}}=8\) for the three examples.

Fig. 1.
figure 1

Distribution of fractures: Example (1), Example (2), Example (3).

Fig. 2.
figure 2

Distribution of the contours and profiles of solute concentration at different dimensionless times 25, 35, 45 and 85.

Fig. 3.
figure 3

Distribution of the contours and profiles of solute concentration at different dimensionless times 15, 35, 55 and 85.

Fig. 4.
figure 4

Distribution of the contours and profiles of solute concentration at different dimensionless times 15, 25, 35 and 60.

The distribution of the contours and profiles of solute concentration in the fractured medium of Example (1) at different dimensionless times 25, 35, 45 and 85 are shown in Fig. 2. This figure indicates that the solute-water mixture moves rapidly in the horizontal fractures due to their high permeability compare to the matrix permeability. The concentration profiles (shown in the right section of Fig. 2) are plotted at three vertical sections at \(y=3.5, 7\) and 10. Again, one may notice that the concentration in fractures is much higher than it in matrix blocks. Similarly, Fig. 3 shows the contours and profiles of solute concentration in the fractured medium of Example (2) at different dimensionless times 15, 35, 55 and 85. One may observe that the mixture flow in fractures is higher than it in the matrix blocks. Finally, the distribution of the contours and profiles of solute concentration at different dimensionless times 15, 25, 35 and 60 of Example (3) are presented in Fig. 4. Similar interpretation to Examples (1) and (2) can be given for Example (3).

5 Conclusions

In this work, we investigate the solute transport using a multi-scale time-stepping method for a single-phase flow in fractured porous media. For the spatial discretion, we used the CCFD method. In the matrix domain a large time-step is used, while in fractures a smaller one is used. The time-step of the fracture pressure is divided into smaller sub-steps, so we explicitly update the concentration at the same time-discretization level. Likewise, the concentration has a bigger time-step in the matrix blocks, while finer ones in the fractures. Three numerical examples were provided to show the efficiency of the proposed scheme. We found that the water-solute quickly transfers through fractures. Therefore, existing fractures in the aquifer improves the distribution of the solute.