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., , and , respectively. In a heterogeneous matrix, the diffusion equation can be written as:
The grey level of each voxel of the image can be converted into an effective porosity value and the tortuosity of diffusion path is assumed to be expressed by an inverse power of the immobile porosity . The boundary condition at the domain boundary between mobile/immobile zone and initial condition . The spatially variable diffusion coefficient in the immobile region depends on the microporosity and tortuosity such that
The tortuosity 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 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.
. The value of is known, so we assume the value of is up or equal to . Subsequently the diffusion in the immobile domain is given simply given by
where 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
with the probability corresponding to the fraction of the volumetric flux that passes by diffusion from i to j such as:
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 by the porosity 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:
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 up to 108 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 . One of three possible movements is executed: 1) The walker is free to move. 2) the walker strikes a nonpermeable solid grain ( bigger than a given threshold ). 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 of the immobile domain. Until the last particle leaves the immobile domain, the total number of random walkers inside the immobile domain is recorded as function of time. Finally, we determine the memory function by the formula given by (Gouze et al., 2008)
with being the diffusion scale defined by
And is a small number with 0 < a < 1 and L the length of a voxel.
4. Simulation Results and Discussion
4.1. Memory Function GXMT
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 , and the tortuosity in the immobile domain. The memory function has been computed for all sections and 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, at early time . 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 (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
in 3D heterogeneous matrix structure using the XMT images (microscale approach). Here, we quantify the roles of the porosity at the percolation threshold
on the behaviour of memory function curves. Results are compared (Figure 3) for different
values (0, 0.4, 0.5 and 0.6) and show that the memory function slopes
and the late-time breakthrough are not similar. At early time, the memory function slopes decrease when
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
the slope of the memory function is very close to the value measured on the tracer test BTCs at intermediate time
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 and 0.6A power law relationship 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 GXMT Versus GMIM
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 C0, the concentration of the injected tracer, and the injection duration, . The unusual shape at the end of the
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 GMIM 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 GMIM is compared to GXMT computed from a section. The two memory functions are similar at the beginning ( , which emphasizes that the scaling of the memory function GMIM is well explained by pore scale diffusion processes for the shortest diffusion times. However, the value of the transition time between behavior and behavior is noticeably smaller (almost one decade) for the pore-scale computed memory function GXMT. 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 GMIM obtained from the SWIW tracer tests with that computed in the subvolume V and cross section S (Figure 5). The computed slope for GXMT at
and n = 6 in 2D and 3D is similar to the value measured on the tracer test BTCs at an early time
for t < 103. However, for 103
Figure 5. Memory function GMIM from the best fit of the SWIW tracer test BTC compared to the memory function GXMT computed using the XMT cross sub-volume (3D) and cross section (2D) in sample S94_1_1.
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 and 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.
The authors are thankful to Lebanese University for the financial support.
BTCs: Break Through Curve
XMT: X-Ray Micro Tomography
TDRW: Time Diffusion Random Walk
SWIW: Single Well Injection Withdrawl
CTRW: Continuous Time Random Walk
 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.
 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.
 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.
 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.
 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.
 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.
 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.
 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.
 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.
 Levy, M., & Berkowitz, B. (2003). Measurement and Analysis of Non-Fickian Dispersion in Heterogeneous Porous Media. Journal of Contaminant Hydrology, 64, 203-226.
 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.
 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.
 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.
 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.
 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.
 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.
 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.