^{1}

^{1}

^{1}

^{1}

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.

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.

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 (

In the immobiledomain, the complicated microstructure causes the microporosity, tortuosity, and the diffusion coefficient, to vary spatially i.e., ϕ ′ i m = ϕ ′ i m ( x ′ ) , τ i m = τ i m ( x ′ ) and d i m = d i m ( x ′ ) , respectively. In a heterogeneous matrix, the diffusion equation can be written as:

ϕ ′ i m ∂ c i m ( x ′ , t ) ∂ t = ∇ ⋅ d i m ( x ′ ) ∇ c i m ( x ′ , t ) (1)

The grey level of each voxel of the image can be converted into an effective porosity value ϕ ′ i m ( x ′ ) and the tortuosity of diffusion path is assumed to be expressed by an inverse power of the immobile porosity τ ( x ′ ) = ϕ ′ ( x ′ ) − n . The boundary condition c i m ( x ′ , t ) | x ′ ϵ ∂ Ω = C m ( x , t ) at the domain boundary ∂ Ω between mobile/immobile zone and initial condition c i m ( x ′ , t = 0 ) | = 0 . The spatially variable diffusion coefficient d ( x ′ ) in the immobile region depends on the microporosity ϕ ′ ( x ′ ) and tortuosity τ ( x ′ ) such that

d i m ( x ′ ) = d 0 ϕ ′ i m ( x ′ ) / τ i m ( x ′ ) (2)

The tortuosity τ 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 τ i m of diffusion path is assumed to be expressed by an inverse power of the immobile porosity

τ i m ( x ′ ) = ϕ ′ i m ( x ′ ) − n . The value of < ϕ ′ i m ( x ′ ) > is known, so we assume the value of τ i m is up or equal to τ m . Subsequently the diffusion in the immobile domain is given simply given by

d i m ( x ′ ) = d 0 ϕ ′ i m ( x ′ ) 1 + n (3)

where d 0 is the molecular diffusion coefficient.

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

Γ i → j = − V i ∑ j b i j (4)

with the probability W i j corresponding to the fraction of the volumetric flux that passes by diffusion from i to j such as:

W i → j = b i j ∑ j b i j (5)

with b i j = S i j ( d i m ) i j / L i j .

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.

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 ∂ C m / ∂ t by the porosity ϕ i m 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 i m = ϕ i m ∂ C m ∂ t ∗ G ( x , t ) (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^{8} 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 Ω = ε S Ω . One of three possible movements is executed: 1) The walker is free to move. 2) the walker strikes a nonpermeable solid grain ( ϕ ′ i m 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 S Ω of the immobile domain. Until the last particle leaves the immobile domain, the total number of random walkers inside the immobile domain N V ( t ′ ) is recorded as function of time. Finally, we determine the memory function by the formula given by (Gouze et al., 2008)

G ( t ) = S Ω V Ω ε N V ( t τ d ) τ d N s (7)

with τ d being the diffusion scale defined by

τ d = ε 2 / ( 2 D d 0 ) (8)

And ε = a L is a small number with 0 < a < 1 and L the length of a voxel.

First, the aim is to compare the memory function in section and volume from the sample presented in

In this section, we present the simulation results of the memory function G X M T 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 (

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.

_{0}, the concentration of the injected tracer, and the injection duration, Δ t . 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). _{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 ( ( t < 10 3 ) , 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 ( t ) ∼ t − 1 behavior and G ( t ) ∼ 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 (_{XMT} at ξ ϕ = 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 ( t ) ∼ t − 1 for t < 10^{3}. However, for 10^{3} 4 the G ( t ) ∼ t − 0.5 . At a late time t > 10 ^{4} the adjusted value of n = 6 corresponds to a Higher tortuosity τ i m 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 τ m in the mobile domain. Accordingly, when τ i m > τ 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 i m ∼ ϕ n .

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 τ i m ∼ ϕ i m − n and ξ ϕ = 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.

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

The authors declare no conflicts of interest regarding the publication of this paper.

Abou-Saleh, K., Dweik, J., Haidar, Y., & Ghaddar, A. (2019). Solving Diffusion Time in Heterogeneous Microscale Rock Matrix by 3D Computations: Non-Fickian Dispersion Observed in Porous Media. Journal of Geoscience and Environment Protection, 7, 42-52. https://doi.org/10.4236/gep.2019.712003

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