Solving Diffusion Time in Heterogeneous Microscale Rock Matrix by 3D Computations: Non-Fickian Dispersion Observed in Porous Media

Khalil Abou-Saleh^*, Jalal Dweik, Youssef Haidar, Abbas Ghaddar

Department of Physics, Faculty of Sciences V, Lebanese University, Nabatieh, Lebanon.

DOI: 10.4236/gep.2019.712003                

Abstract

Understanding and simulating the underlying microscopic physics of the rock matrix is very useful for determining macroscopic physical properties such as permeability. Matrix diffusion is an important transport parameter controlling the late-time behaviour of breakthrough curves (BTCs). We compute the memory function, implemented in the sink/source term of Mobile-immobile mass transfer by solving the matrix diffusion using a time diffusion random-walk approach. The diffusion is controlled by different parameters like the porosity, tortuosity, mobile-immobile interface and immobile domain cluster shapes. All these properties are investigated by X-ray microtomography that captures the main characteristics of matrix diffusion at three dimensions. We compare the memory function deduced from the field-scale tracer tests well with the computed memory function. Simulation results of the memory function appeared to be coherent with that measured from the tracer test for a large tortuosity value. Probably, the diffusion paths are longer, and they are controlled by the properties mentioned above. From a representative elementary volume of natural reservoirs studied here, we conclude that, microscale diffusion process in the immobile domain play a crucial role to better understand the non-Fickian dispersion measured from the tracer test.

Keywords

Matrix Diffusion, Time Domain Random Walk, Memory Function, Mobile-Immobile Mass Transfer

1. Introduction

Transport properties in porous media are probably the most important parameters in many geophysical and engineering situations, including pollutant migration. Porous media, such as rocks, consist of pore space and a solid matrix. The pore spaces are usually connected, which allows fluid flow and mass transfer to take place. Permeability and electrical conductivity, as well as other parameters, are strongly influenced by the pore structure and pore-scale physics. In fact, we can determine the exact transport properties (Cai et al., 2019) of a porous medium if we can solve the physical equations on a real pore space. However, Matrix diffusion has an effective effect on the transport mechanism (Zhan et al., 2019; Zhou et al., 2017), controlling the late-time behaviour of breakthrough curves (BTCs) (Zhou et al., 2018; Shapiro, 2001; Zhoua et al., 2007; Gouze et al., 2008). An immobile volume zone in the matrix becomes accessible to solutes by diffusion. This process causes a delay of solutes arrival times and consequently produces strongly asymmetric BTCs (non-Fickian) at the macroscopic scale (Zhoua et al., 2007; Gouze et al., 2008; Berkowitz & Scher, 2000; Meigs & Beauhelm, 2001; Levy & Berkowitz, 2003). Many authors focus on the transport in the rock’s matrix. Carrera (Carrera et al., 1998) proposed a formulation of the sink/source term of the transport equation, to characterize diffusion in the simple geometry of matrix (e.g., labs and spheres). Reference (Haggerty et al., 2000; Haggerty et al., 2001; Haggerty et al., 2004) reported a study on the transport in the rock matrix using the multiple rates of mass transfers for any complex structures (Gouze et al., 2008) reported a study to solve the transport using a 2-D microscale description of the matrix structure, precisely by computing the memory function in Mallorc limestones sample using the classical random walk particle tracking method (Salamon et al., 2006; Delay et al., 2005). These authors were inspired by the non-Fickian dispersion observed on tracer test BTCs. The idea is to link this observation from field scale, to microscale diffusion process within water immobile zones. The results obtained by (Gouze et al., 2008) showed that large-scale non-Fickian dispersion is controlled by the microscale matrix diffusion.

This class of model is often referred to as a mobile/immobile (MIM) mass transfer model by which solutes transfer from the water-flowing portions (mobile) of permeable media to the non-flowing portions (immobile). In recent years, computing power is increasing; the simulation techniques are becoming an alternative solution. Recently, (Dweik et al., 2015) showed that computed memory function from three-dimensional XMT images provides a more accurate definition than the two-dimensional ones. The problem that should be addressed is the delayed time behavior of BTCs that is related to the memory. The motivation of this study is the results obtained (Gouze et al., 2008). From XMT cross sections, the residence time of the tracer in the immobile domain is underestimated when computed, may be the calculations carried out on longer cross sections diffusion paths. Probably three-dimensional computations tend to lengthen the time axis if the diffusion paths are longer.

In this paper, we determine the memory function by 3D computations in microscale matrix structure, whose properties are investigated by X-ray microtomography (XMT). We solve the transport in heterogeneous matrix by the time diffusion random walk (TDRW) method inspired by (Delay et al., 2002; Sardini et al., 2003; Sardini et al., 2007). Then, we develop an algorithm to compute the memory function in real matrix structure extracted from X-ray micro-tomo- graphy (XMT). In addition, the tortuosity parameter and its effect on the behaviour of memory function at the different percolation thresholds is quantified. Finally, we compare the computed memory function with those deduced from the field-scale tracer tests.

2. Experimental Data

The tracer tests discussed here are performed in a Miocene reef formation situated 50 km from Palma de Mallorca (Balearic Islands). The reservoir rock is a pure bioclastic carbonate (calcite). The characteristics of the medium, the tracer test equipment and protocols as well as the results obtained for the 10 tests performed at depth, 94 m are detailed by (Gouze et al., 2008; Khrapitchev & Callaghan, 2003). Cylindrical minicores of 18 mm length and 10 mm diameter were sampled from the borehole core at a depth corresponding to the SWIW tracer tests. Micro-Tomography XMT was used to investigate the pore space of Mallorca limestone. Data were collected using the BM5 beam line of the European Synchrotron Radiation Facility (Grenoble, France). The stages of identification of different regions in the cluster (Figure 1) are well detailed in the previous article (Dweik et al., 2015).

3. Mathematical Formulations and Algorithm

In the immobiledomain, the complicated microstructure causes the microporosity, tortuosity, and the diffusion coefficient, to vary spatially i.e., ${{\varphi }^{\prime }}_{im}={{\varphi }^{\prime }}_{im}\left(x^{\prime }\right)$, ${\tau }_{im}={\tau }_{im}\left(x^{\prime }\right)$ and $d_{im}=d_{im}\left(x^{\prime }\right)$, respectively. In a heterogeneous matrix, the diffusion equation can be written as:

${{\varphi }^{\prime }}_{im}\frac{\partial c_{im}\left(x^{\prime },t\right)}{\partial t}=\nabla \cdot d_{im}\left(x^{\prime }\right)\nabla c_{im}\left(x^{\prime },t\right)$ (1)

The grey level of each voxel of the image can be converted into an effective porosity value ${{\varphi }^{\prime }}_{im}\left(x^{\prime }\right)$ and the tortuosity of diffusion path is assumed to be expressed by an inverse power of the immobile porosity $\tau \left(x^{\prime }\right)={\varphi }^{\prime }{\left(x^{\prime }\right)}^{-n}$. The boundary condition ${c_{im}\left(x^{\prime },t\right)|}_{x^{\prime }ϵ\partial \Omega }=C_{m\left(x,t\right)}$ at the domain boundary $\partial \Omega$ between mobile/immobile zone and initial condition $c_{im}\left(x^{\prime },t=0\right)|=0$. The spatially variable diffusion coefficient $d\left(x^{\prime }\right)$ in the immobile region depends on the microporosity ${\varphi }^{\prime }\left(x^{\prime }\right)$ and tortuosity $\tau \left(x^{\prime }\right)$ such that

$d_{im}\left(x^{\prime }\right)=d_0{{\varphi }^{\prime }}_{im}\left(x^{\prime }\right)/{\tau }_{im}\left(x^{\prime }\right)$ (2)

The tortuosity ${\tau }_m$ is an important transport parameter controlling the arrival time of the tracers in porous media (mobile zone). These parameters can be deduced from numerical simulations of time-dependent effective diffusion coefficient in the porosity of the rock. In the matrix, the tortuosity ${\tau }_{im}$ of diffusion path is assumed to be expressed by an inverse power of the immobile porosity

Figure 1. Three-dimensional volume (510 × 510 × 510 pixels) that represents of a sub-volume in sample with the coordinates are given in pixels (pixel size 5.01). The immobile fraction is displayed in blue, the mobile fraction of the connected macroporosity is displayed in cyan, the porosity of the microporous phase is displayed in shades of gray with dark and light gray denote low and high porosity, respectively, and the non-connected macroporosity is displayed in green.

${\tau }_{im}\left(x^{\prime }\right)={{\varphi }^{\prime }}_{im}{\left(x^{\prime }\right)}^{-n}$. The value of $<{{\varphi }^{\prime }}_{im}\left(x^{\prime }\right)>$ is known, so we assume the value of ${\tau }_{im}$ is up or equal to ${\tau }_m$. Subsequently the diffusion in the immobile domain is given simply given by

$d_{im}\left(x^{\prime }\right)=d_0{{\varphi }^{\prime }}_{im}{\left(x^{\prime }\right)}^{1+n}$ (3)

where $d_0$ is the molecular diffusion coefficient.

3.1. Time Domain Random Walk (TDRW)

The TDRW simulates solute transport in 3D heterogeneous media and calculates the arrival time of a particle cloud at a given location. The principle is to calculate the diffusion time needed by a particle to jump from one center to the other center of neighbouring voxels with the transition probability P. Each particle holds its location and residence time, which are recorded at each jump. This method is inspired by and modified from attempts in petroleum engineering to solve the steady-state flow equation (Delay et al., 2002; McCarthy, 1993; Dentz et al., 2012). The diffusion time from i to j across the volume of voxel I is given by

${\Gamma }_{i\to j}=-\frac{V_i}{{\displaystyle {\sum }_j{b}_{ij}}}$ (4)

with the probability ${\mathcal{W}}_{ij}$ corresponding to the fraction of the volumetric flux that passes by diffusion from i to j such as:

${\mathcal{W}}_{i\to j}=\frac{b_{ij}}{{\displaystyle {\sum }_j{b}_{ij}}}$ (5)

with $b_{ij}=S_{ij}{\left(d_{im}\right)}_{ij}/L_{ij}.$

The Equations (4) and (5) can be integrated in the algorithm of memory function to determine respectively the diffusion time and the displacement of the particles.

3.2. Algorithm of Memory Function

In the mobile domain, the sink/source term implemented by the transport equation can be expressed as the convolution product of the time variation of the mobile concentration $\partial C_m/\partial t$ by the porosity ${\varphi }_{im}$ in immobile domain and time-dependent function, called the memory function G(x,t), (Gouze et al., 2008; Carrera et al., 1998; Haggerty et al., 2000) such that:

$F_{im}={\varphi }_{im}\frac{\partial C_m}{\partial t}\ast G\left(x,t\right)$ (6)

the memory function is an intrinsic characteristic of the medium, depends only on the properties of the immobile domain, then the Equation (6) is easy to compute.

The memory function G(t) is the probability that a tracer entering the immobile zone will stay there until time t (Haggerty et al., 2004). We develop an algorithm to compute the memory function using the TDRW method. The goal is to calculate the diffusion time required for a particle to go from one center to another of the neighboring pixels according to the transition probability P. Large number of random walkers $N_s$ up to 10⁸ is applied in the simulation to obtain clear results. The initial distribution of the random walkers at time t = 0, is at the mobile-immobile interface over a small volume $V_{S_{\Omega }}=\epsilon S_{\Omega }$. One of three possible movements is executed: 1) The walker is free to move. 2) the walker strikes a nonpermeable solid grain ( ${{\varphi }^{\prime }}_{im}$ bigger than a given threshold ${\xi }_{\varphi }$ ). In this case the jump is not performed, and a new random position is selected for a new position to be observed. 3) The walker crosses the surface $S_{\Omega }$ of the immobile domain. Until the last particle leaves the immobile domain, the total number of random walkers inside the immobile domain $N_V\left(t^{\prime }\right)$ is recorded as function of time. Finally, we determine the memory function by the formula given by (Gouze et al., 2008)

$G\left(t\right)=\frac{S_{\Omega }}{V_{\Omega }}\frac{\epsilon N_V\left(\frac{t}{{\tau }_d}\right)}{{\tau }_d{N}_s}$ (7)

with ${\tau }_d$ being the diffusion scale defined by

${\tau }_d={\epsilon }^2/\left(2D{d}_0\right)$ (8)

And $\epsilon =aL$ is a small number with 0 < a < 1 and L the length of a voxel.

4. Simulation Results and Discussion

4.1. Memory Function G_XMT

First, the aim is to compare the memory function in section and volume from the sample presented in Figure 1. The residence time of the solute in the immobile domain is computed, using the algorithm described in Section 3. All the parameters or the properties extracted from the samples (sub-volume V) are known except for the value of the porosity at the percolation ${\xi }_{\varphi }$, and the tortuosity in the immobile domain. The memory function has been computed for all sections $S_{xy},S_{yz}$ and $S_{xz}$ in the sub-volume sample. Results are compared with the memory function computed from the total sub-volume V. Figure 2 shows that the memory function G(t) slopes are similar, $\gamma \approx 1$ at early time $t<{10}^2$. However, longer diffusion time in three-dimensional computations is observed than in two-dimensions at late time. This might be linked to several factors. First, the particles in 3D walk a long path controlled by the effective value of the porosity at the percolation threshold. In terms of a continuous model, forming a connected voxel belonging to the immobile domain is equal to its cluster of pores and requires porosity above a given threshold ${\xi }_{\varphi }$ (Gouze, 2008). Also, the matrix block shape is integrated with the surface. Furthermore, the heterogeneity of the porosity is in different directions (x, y, z). 3D computations stretch the time axis because the diffusion paths are longer. The memory function computed from XMT images at three-dimensions provides a more accurate definition rather than the two-dimensional, as the simulations results, shows a longer residence time of particles in the immobile domain.

In this section, we present the simulation results of the memory function $G_{XMT}$ in 3D heterogeneous matrix structure using the XMT images (microscale approach). Here, we quantify the roles of the porosity at the percolation threshold ${\xi }_{\varphi }$ on the behaviour of memory function curves. Results are compared (Figure 3) for different ${\xi }_{\varphi }$ values (0, 0.4, 0.5 and 0.6) and show that the memory function slopes $\gamma$ and the late-time breakthrough are not similar. At early time, the memory function slopes decrease when ${\xi }_{\varphi }$ increases. The diffusion path of the particles is controlled by the surface block shape (mobile-immobile interface geometry) of the matrix and by the heterogeneity of the porosity. For ${\xi }_{\varphi }=0.6$ the slope of the memory function is very close to the value measured on the tracer test BTCs at intermediate time $t_{ad}<t<t_c$ ( $t_{ad}$ and are the advection and transition time respectively). In contrast, for a long time, the residence time of the particles becomes longer when ${\xi }_{\varphi }$ increases, according to the probability of displacements. In this context, we quantify the roles of the tortuosity $\tau$. This factor is certainly the most important parameter controlling the effective diffusion of the tracers without forgetting the mobile-immobile interface geometry. A power law relationship, ${\tau }_{im}~{\varphi }^{-n}$ with different value of n = 0, 2 and 4, is used in the numerical simulation. The influence of the tortuosity

Figure 2. Computed memory function G(t) in 2-D structure with 510 × 510 pixels size along z direction versus G(t) in 3-D structure with 510 × 510 × 510 pixels size.

Figure 3. Memory function computed in cross sub-volume V for different value of ${\zeta }_{\varphi }=0,0.4,0.5$ and 0.6A power law relationship ${\tau }_{im}\sim {\varphi }^{-n}$ at different value of n = 0, 2 and 4 is used in the numerical simulation.

value on the memory function shape appeared on the late time breakthrough. It is clear that the residence time of the particles inside the matrix increase when the value of tortuosity increases.

4.2. Comparison of G_XMT Versus G_MIM

Figure 4(a) shows the Breakdown Curve (BTC) obtained from the plotting experiment and the fitted curves of the MIM model (Gouze et al., 2008). Concentration, C, is normalized by dividing it by C₀, the concentration of the injected tracer, and the injection duration, $\Delta t$. The unusual shape at the end of the

Figure 4. (a) Best fit of the SWIW tracer test BTC; (b) Memory function G_XMT computed using the XMT cross section S4 compared to the memory function (G_MIM) (Gouze et al., 2008).

curve is modeled by an MIM mass transfer model using the CTRW approach with a double slope transition time distribution (Le Borgne & Gouze, 2008). Figure 4 (right), shows the G_MIM corresponding to the SWIW tracer tests, such as the memory function best fits of BTC, is obtained by a trial-and-error method (Carrera et al., 1998). The G_MIM is compared to G_XMT computed from a section. The two memory functions are similar at the beginning ( $\left(t<{10}^3\right)$, which emphasizes that the scaling of the memory function G_MIM is well explained by pore scale diffusion processes for the shortest diffusion times. However, the value of the transition time between $G\left(t\right)\sim t^{-1}$ behavior and $G\left(t\right)\sim t^{-0.5}$ behavior is noticeably smaller (almost one decade) for the pore-scale computed memory function G_XMT. Similarly, the cutoff time is about two decades less than the minimum cutoff time required to fit the tracer test data and the effective value of the porosity at the percolation threshold should be different (Gouze et al., 2008).

In this section, we compare the memory function G_MIM obtained from the SWIW tracer tests with that computed in the subvolume V and cross section S (Figure 5). The computed slope for G_XMT at ${\xi }_{\varphi }=0.6$ and n = 6 in 2D and 3D is similar to the value measured on the tracer test BTCs at an early time $G\left(t\right)\sim t^{-1}$ for t < 10³. However, for 10³ 4 the $G\left(t\right)\sim t^{-0.5}$. At a late time t > 10 ⁴ the adjusted value of n = 6 corresponds to a Higher tortuosity ${\tau }_{im}$ in the matrix, the memory G _XMT computed from 3D subvolume provides a more accurate definition as mention before. Considering that both the dual slope behaviour is controlled by both the matrix block shape and the heterogeneity of the porosity control the dual slope behaviour, 3D computations stretch the time axis because the diffusion paths are longer. The value of n can be estimated, based on the tortuosity ${\tau }_m$ in the mobile domain. Accordingly, when ${\tau }_{im}>{\tau }_m$ the best value of n that fits the tracers test is equal to 6. For instance, the tortuosity of the diffusion paths, as well as the variability and the spatial correlation of the diffusivity along the paths are probably imperfectly modelled by a relation of the form $d_{im}\sim {\varphi }^n$.

Figure 5. Memory function G_MIM from the best fit of the SWIW tracer test BTC compared to the memory function G_XMT computed using the XMT cross sub-volume (3D) and cross section (2D) in sample S94_1_1.

5. Conclusion

In this paper, we determine the memory function to characterize the diffusion in the matrix (i.e., the immobile domain). This study gives a deep understanding on the large-scale non-Fickian dispersion that is, controlled by the microscale properties of the matrix. The XMT technique appears to be a promising tool to capture the properties of the immobile domain at 3D and then compute the memory function using the time domain random walk. The variation of the memory function at late time is controlled by the tortuosity of the diffusion path and the value of the porosity at the percolation threshold. Using the empirical relation ${\tau }_{im}\sim {\varphi }_{im}^{-n}$ and ${\xi }_{\varphi }=0.6$ for n = 6 the memory function computed from XMT images at three-dimensional provides a more accurate definition rather than the two-dimensional, as the simulations results, show a longer residence time of particles in the immobile domain. Using 3D computation of the memory function can provide insights on the delay of the tracer’s arrival time caused by the matrix structure parameter. Furthermore, the simulation results obtained in this study is coherent with those measured from the tracer test. A diffusion mobile-immobile mass transfer model appears to be valid for the reservoir rock studied, and the atypical non-Fickian dispersion measured from the tracer test well explained by the microscale diffusion processes in the immobile domain. Finally, it is interesting to determine the memory function in another type of reservoir rock.

Acknowledgements

The authors are thankful to Lebanese University for the financial support.

Symbols/Acronyms List

BTCs: Break Through Curve

MIM: Mobile/ImMobile

XMT: X-Ray Micro Tomography

TDRW: Time Diffusion Random Walk

SWIW: Single Well Injection Withdrawl

CTRW: Continuous Time Random Walk

Conflicts of Interest

The authors declare no conflicts of interest.

References

[1]	Berkowitz, B., & Scher, H. (2000). Anomalous Transport in Laboratory-Scale, Heterogeneous Porous Media. Water Resources Research, 36, 149-158. https://doi.org/10.1029/1999WR900295
[2]	Cai, J., Zhang, Z., Kang, Q. et al. (2019). Recent Advances in Flow and Transport Properties of Unconventional Reservoirs. Energies, 12, 1865-1869. https://doi.org/10.3390/en12101865
[3]	Carrera, J., Sánchez-Vila, X., Benet, I. et al. (1998). On Matrix Diffusion, Formulations, Solutions Methods and Qualitative Effects. Hydrogeology Journal, 6, 178-190. https://doi.org/10.1007/s100400050143
[4]	Delay, F., Ackerer, P., & Danquigny, C. (2005). Simulating Solute Transport in Porous or Fractured Formations Using Random Walk Particle Tracking: A Review. Vadose Zone Journal, 4, 360-379. https://doi.org/10.2136/vzj2004.0125
[5]	Delay, F., Porel, G., & Sardini, P. (2002). Modelling Diffusion in a Heterogeneous Rock Matrix with a Time-Domain Lagrangian Method and an Inversion Procedure. C.R. Geoscience, 334, 967-973. https://doi.org/10.1016/S1631-0713(02)01835-7
[6]	Dentz, M., Gouze, P., Russian, A. et al. (2012). Diffusion and Trapping in Heterogeneous Media: An Inhomogeneous Continuous Time Random Walk Approach. Advances in Water Resources, 49, 13-22. https://doi.org/10.1016/j.advwatres.2012.07.015
[7]	Dweik, J., Srour, M., Tawbe, M. et al. (2015). Solving Microscale Heterogeneous Matrix Diffusion Based on Two and Three-Dimensional Computing Using X-Ray Tomography Image Data. International Review of Physics, 9, 79-84.
[8]	Gouze, P., Melean, Y., Le Borgne, T. et al. (2008). Non-Fickian Dispersion in Porous Media Explained by Heterogeneous Microscale Matrix Diffusion. Water Resources Research, 44, W11416-W11435. https://doi.org/10.1029/2007WR006690
[9]	Haggerty, R., Fleming, S. W., Meigs, L. C. et al. (2001). Tracer Tests in a Fractured Dolomite: 2. Analysis of Mass Transfer in Single-Well Injection-Withdrawal Tests. Water Resources Research, 37, 1129-1142. https://doi.org/10.1029/2000WR900334
[10]	Haggerty, R., Harvey, C. F., Schwerin, C. F. et al. (2004). What Controls the Apparent Timescale of Solute Mass Transfer in Aquifers and Soils? A Comparison of Experimental Results. Water Resources Research, 40, W01510. https://doi.org/10.1029/2002WR001716
[11]	Haggerty, R., McKenna, S. A., & Meigs, L. C. (2000). On the Late-Time Behavior of Tracer Test Breakthrough Curves. Water Resources Research, 36, 3467-3480. https://doi.org/10.1029/2000WR900214
[12]	Khrapitchev, A. A., & Callaghan, P.-T. (2003). Reversible and Irreversible Dispersion in a Porous Medium. Physics of Fluids, 15, 2649-2660. https://doi.org/10.1063/1.1596914
[13]	Le Borgne, T., & Gouze, P. (2008). Non-Fickian Dispersion in PorousMedia: 2. Model Validation from Measurements at Different Scales. Water Resources Research, 44, W06427- W06436. https://doi.org/10.1029/2007WR006279
[14]	Levy, M., & Berkowitz, B. (2003). Measurement and Analysis of Non-Fickian Dispersion in Heterogeneous Porous Media. Journal of Contaminant Hydrology, 64, 203-226. https://doi.org/10.1016/S0169-7722(02)00204-8
[15]	McCarthy, J. F. (1993). Continuous Time Random Walks on Random Media. Journal of Physics A: Mathematical and General, 26, 2495-2503. https://doi.org/10.1088/0305-4470/26/11/004
[16]	Meigs, L. C., & Beauhelm, R. L. (2001). Tracer Tests in a Fractured Dolomite: 1, Experimental Design and Observed Tracer Recoveries. Water Resources Research, 37, 1113- 1128. https://doi.org/10.1029/2000WR900335
[17]	Salamon, P., Fernàndez-Garcia, D., & Gómez-Hernández, J. J. (2006). A Review and Numerical Assessment of the Random Walk Particle Tracking Method. Journal of Contaminant Hydrology, 87, 277-305. https://doi.org/10.1016/j.jconhyd.2006.05.005
[18]	Sardini, P., Delay, F., & Hellmuth, K.-H. (2003). Interpretation of Out-Diffusion Experiments on Crystalline Rocks Using Random Walk Modelling. Journal of Contaminant Hydrology, 61, 339-350. https://doi.org/10.1016/S0169-7722(02)00124-9
[19]	Sardini, P., Robinet, J.-C., & Siitari-Kappi, M. (2007). Direct Simulation of Heterogeneous Diffusion and Inversion Procedure Applied to an Out-Diffusion Experiment. Test Case of Palmottugranite. Journal of Contaminant Hydrology, 93, 21-37. https://doi.org/10.1016/j.jconhyd.2007.01.011
[20]	Shapiro, A. M. (2001). Effective Matrix Diffusion in Kilometer-Scale Transport in Fractured Crystalline Rock. Water Resources Research, 37, 507-522. https://doi.org/10.1029/2000WR900301
[21]	Zhan, H., Wang, Q., & Wen, Z. (2019). Advances in Groundwater Flow and Solute Transport: Pushing the Hidden Boundary. Water, 11, 457-460. https://doi.org/10.3390/w11030457
[22]	Zhou, L., Jing, L., & Cvetkovic, V. (2017). Modeling of Solute Transport in a 3D Rough- Walled Fracture-Matrix System. Transport in Porous Media, 116, 1005-1029. https://doi.org/10.1007/s11242-016-0810-z
[23]	Zhou, Q., Oldenburg, C.-M., Kneafsey, T.-J. et al. (2018). Modeling Transport of Multiple Tracers in Hydraulic Fractures at the, EGS Collab Test Site. In 43rd Workshop on Geothermal Reservoir Engineering (SGP-TR-213). Stanford: Stanford University.
[24]	Zhoua, Q., Liua, H.-H., Molz, F. J. et al. (2007). Field-Scale Effective Matrix Diffusion Coefficient for Fractured Rock: Results from Literature Survey. Journal of Contaminant Hydrology, 93, 161-187. https://doi.org/10.1016/j.jconhyd.2007.02.002