Abstract
In this paper, we investigate flow with solute transport in fractured porous media. The system of the governing equations consists of the continuity equation, Darcy’s law, and concentration equation. A discrete-fracture model (DFM) has been developed to describe the problem under consideration. The multiscale time-splitting method was used to handle different sizes of time-step for different physics, such as pressure and concentration. Some numerical examples are presented to show the efficiency of the multi-scale time-splitting approach.
You have full access to this open access chapter, Download conference paper PDF
Similar content being viewed by others
Keywords
- Multiscale time-splitting
- Discrete fracture model
- Mass transfer
- Fracture porous media
- Reservoir simulation
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):
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:
where q is the external mass flow rate.
Mass Transfer: The solute transport equation in porous media may be given as,
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,
The boundary conditions are given as,
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,
Assuming that the pressure along the fracture width are constants, and by integration, the pressure equation in the fracture becomes,
The matrix-fracture interface condition is given by,
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,
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,
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,
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,
and
Now, applying the CCFD scheme on (13), one obtains,
If \(b\in \partial K\cap \partial K'\) and \(b\nsubseteq \varOmega _{\text {f}}\), the fluxes in (15) are given by,
where \(\chi _{t,b}^i\) is given by
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,
where \(\chi _{t,\text {mf},b}^i\) is given by,
Similarly, let b be a gridcell of the fracture network. Equation (14) may be discretized by the CCFD method to get,
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.
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,
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,
The pressure equation in the fractures at each time subtime-step is given by,
It is obtained by using the CCFD scheme to (26) that
The matrix-fracture transfer is given by,
where
along with,
For the case of multiple fractures, the pressure equation may be given as,
where
and
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}}\),
where
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,
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,
We use the upwind CCFD method to discretize the concentration Eq. (38),
where
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,
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,
where \(\chi _{t,b}^{i,j,k}\) is given by the harmonic mean as,
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,
where \(\chi _{t,\text {mf},b}^{i,j,k,l+1}\) is defined as,
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.
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.
References
Park, Y.J., Lee, K.K.: Analytical solutions for solute transfer characteristics at continuous fracture junctions. Water Resour. Res. 35(5), 1531–1537 (1999)
Roubinet, D., de Dreuzy, J.R., Daniel, M.T.: Semi-analytical solutions for solute transport and exchange in fractured porous media. Water Resour. Res. 48(1), W01542 (2012)
Bodin, J., Porel, G., Delay, F.: Simulation of solute transport in discrete fracture networks using the time domain random walk method. Earth Planet. Sci. Lett. 208(3–4), 297–304 (2003)
Graf, T., Therrien, R.: Variable-density groundwater flow and solute transport in irregular 2D fracture networks. Adv. Water Resour. 30(3), 455–468 (2007)
Jackson, C.P., Hoch, A.R., Todman, S.: Self-consistency of a heterogeneous continuum porous medium representation of a fractured medium. Water Resour. Res. 36(1), 189–202 (2000)
Long, J.C.S., Remer, J.S., Wilson, C.R., Witherspoon, P.A.: Porous media equivalents for networks of discontinuous fractures. Water Resour. Res. 18(3), 645–658 (1982)
Weatherill, D., Graf, T., Simmons, C.T., Cook, P.G., Therrien, R., Reynolds, D.A.: Discretizing the fracture-matrix interface to simulate solute transport. Ground Water 46(4), 606–615 (2008)
McKenna, S.A., Walker, D.D., Arnold, B.: Modeling dispersion in three-dimensional heterogeneous fractured media at Yucca Mountain. J. Contam. Hydrol. 62–63, 577–594 (2003)
Zhou, Z.F., Guo, G.X.: Numerical modeling for positive and inverse problem of 3-D seepage in double fractured media. J. Hydrodyn. 17(2), 186–193 (2005)
Kou, J., Sun, S., Yu, B.: Multiscale time-splitting strategy for multiscale multiphysics processes of two-phase flow in fractured media. J. App. Math. 2011, Article ID 861905, 24 pages (2011)
El-Amin, M.F., Khalid, U., Beroual, A.: Magnetic field effect on a ferromagnetic fluid flow and heat transfer in a porous cavity. Energies 11(11), 3235 (2018)
El-Amin, M.F., Kou, J., Sun, S., Salama, A.: Adaptive time-splitting scheme for two-phase flow in heterogeneous porous media. Adv. Geo-Energy Res. 1(3), 182–189 (2017)
El-Amin, M.F., Kou, J., Sun, S.: Adaptive time-splitting scheme for nanoparticles transport with two-phase flow in heterogeneous porous media. In: Shi, Y., et al. (eds.) ICCS 2018. LNCS, vol. 10862, pp. 366–378. Springer, Cham (2018). https://doi.org/10.1007/978-3-319-93713-7_30
El-Amin, M.F., Kou, J., Sun, S.: Discrete-fracture-model of multi-scale time-splitting two-phase flow including nanoparticles transport in fractured porous media. J. Comp. App. Math. 333, 327–249 (2018)
Author information
Authors and Affiliations
Corresponding author
Editor information
Editors and Affiliations
Rights and permissions
Copyright information
© 2020 Springer Nature Switzerland AG
About this paper
Cite this paper
El-Amin, M.F., Kou, J., Sun, S. (2020). Numerical Investigation of Solute Transport in Fractured Porous Media Using the Discrete Fracture Model. In: Krzhizhanovskaya, V., et al. Computational Science – ICCS 2020. ICCS 2020. Lecture Notes in Computer Science(), vol 12143. Springer, Cham. https://doi.org/10.1007/978-3-030-50436-6_8
Download citation
DOI: https://doi.org/10.1007/978-3-030-50436-6_8
Published:
Publisher Name: Springer, Cham
Print ISBN: 978-3-030-50435-9
Online ISBN: 978-3-030-50436-6
eBook Packages: Computer ScienceComputer Science (R0)