<figure><img src="image_1.png" /></figure>

HAL
open science

# Multi-Scale Study of Diffusion in Composite Grain-Pore Systems Based on Particles Random Walk

Hamza Oukili, Rachid Ababou, Gerald Debenest, Benoît Noetinger

► To cite this version:

Hamza Oukili, Rachid Ababou, Gerald Debenest, Benoît Noetinger. Multi-Scale Study of Diffusion in Composite Grain-Pore Systems Based on Particles Random Walk. Comptes Rendus. Mécanique, 2021, 349 (3), pp.529-558. 10.5802/crmeca.94. hal-03499827

HAL Id: hal-03499827

https://ifp.hal.science/hal-03499827

Submitted on 21 Dec 2021

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_2.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# Comptes Rendus
## Mécanique

Hamza Oukili, Rachid Ababou, Gérald Debenest and Benoît Noetinger

**Multi-scale study of diffusion in composite grain-pore systems based on particles random walk**

Volume 349, issue 3 (2021), p. 529-558

<https://doi.org/10.5802/crmeca.94>

© Académie des sciences, Paris and the authors, 2021.
*Some rights reserved.*

<figure><img src="image_3.png" /></figure>

This article is licensed under the
CREATIVE COMMONS ATTRIBUTION 4.0 INTERNATIONAL LICENSE.
http://creativecommons.org/licenses/by/4.0/
<figure><img src="image_4.png" /></figure>

Short paper / *Note*

# Multi-scale study of diffusion in composite grain-pore systems based on particles random walk

Hamza Oukili<sup>●</sup> *,<sup>a</sup> Rachid Ababou<sup>a</sup>, Gérald Debenest<sup>a</sup>
and Benoît Noetinger<sup>● b</sup>

<sup>a</sup> Institut de Mécanique des Fluides de Toulouse, IMFT, Université de Toulouse, CNRS
- Toulouse, France

<sup>b</sup> IFP Energies nouvelles, Rueil-Malmaison, France

*E-mails:* hamza.oukili.01@gmail.com (H. Oukili), rachidababou@gmail.com (R. Ababou), gerald.debenest@toulouse-inp.fr (G. Debenest), benoit.noetinger@ifpen.fr (B. Noetinger)

**Abstract.** This paper studies diffusion processes at different scales in two-dimensional (2D) composite media (grain-pore porous micromodels). To model diffusion, this study uses the random walk particle tracking (RWPT) method based on 2D analytical solutions in finite and semi-infinite domains with zero-flux boundary conditions (BC) at grain boundaries. The analytical solutions are developed using a method of images. Then, an RWPT combined to a Stop&Go algorithm is derived from these analytical solutions. The developed RWPT algorithm is then used to model diffusion processes inside the pores at the “microscopic” level (microscale), while the grain elements are assumed inaccessible to diffusion (zero-flux BC condition at the interface grain/pore). The composite medium, where the diffusion occurs, is a numerical micromodel made of periodic motifs of rectangular grains. Particles are initially located at the center of a pore, and diffuses in the periodic motifs micromodel (infinite domain). First, the grains are chosen square with different sizes to vary the porosity. Second, for constant porosities, the grains are elongated to study the effect of the anisotropy ratio on diffusion processes. Effective macroscale properties (porosities, effective diffusion tensors, tortuosities) are then calculated using moments of particles positions. The results obtained fit well with theoretical expectations and are in very good agreement with results found in the literature.

**Keywords.** Composite media, Diffusion, Random walk particle tracking (RWPT), Discontinuities, Analytical solutions, Porous materials, Zero flux.

*Manuscript received 7th February 2021, revised 10th October 2021, accepted 6th September 2021.*

## 1. Introduction

Transport processes occurs in many fields of science including, mass or heat transfer problem, for instance. In porous media, they include also solute mass transport, heat transport, and slightly

<footnote>* Corresponding author.</footnote>
compressible Darcy flow (pressure diffusion equation). Diffusion processes, in particular, are
important in transport phenomena, either directly (solute diffusion, heat diffusion, pressure
diffusion) or indirectly (when the local hydrodynamic dispersion of transported quantities is
modeled as a diffusion-type process on a larger scale). Treating such problems in physics could
be challenging, particularly in the case of spatially variable or discontinuous transport properties.

At the pore scale, considering a binary grain/pore medium, diffusion usually occurs inside the pores, while the grain elements are assumed inaccessible to diffusion (extensions to aggregates instead of grains are discussed in the conclusion Section 5). A zero-flux condition is assumed locally at the grain/pore interfaces. When we treat the diffusion process at the macroscale level, diffusion occurs in a homogeneous medium with macroscopic parameters (including porosity and effective diffusion coefficients) determined via an upscaled methodology. There exist several techniques to determine correctly the link between the transport processes at the microscale and the effective parameters at the macroscale. In the literature, diffusion problems and their upscaling issues have been widely treated using either volume averaging techniques [1, 2], homogenization theory [3, 4], or moments matching techniques [5, 6]. A review on multi-scale methods [7] has shown that homogenization still faces many challenges. Another recent study has shown some improvements by using different levels of upscaling in modeling flow transport in fractured media [8]. Upscaling has also been used to determine parameters for continuous time random walk (CTRW) modeling [9]. In addition, the random walk particle methods have also been used to obtain probability results on self-avoiding Random Walks where the “particles” cannot cross their own path [10].

Solving the diffusion equation at various scales is possible using classical methods such as Finite Volume or Finite Elements methods. In this study, we use particle methods. They have been extensively used for modeling transport problems in porous soils, aquifers, and reservoirs. They reduce or avoid some of the problems of Eulerian methods, for example, instabilities, excessive artificial diffusion, mass balance, and/or oscillations that could lead to negative concentrations. In the Lagrangian particle approach, diffusion processes are modeled by random walk algorithms. The numerical (modeling) issues are different for the Eulerian partial differential equations (PDE) approach and for the Lagrangian random walk approach (for more details, see the recent study of [11] and the review paper of [12]).

Previous studies have used CTRW algorithms for upscaling double porosity fractured media [13–15]. Random walk methods have also been used in many publications in order to analyze the “effective” or “macroscale” diffusion coefficient that emerges from the transport of Brownian particles representing a passive tracer (heat or mass) transported in a heterogeneous or random velocity field (e.g., [16], among others). In this study, we develop a new algorithm, using the random walk particle tracking (RWPT) method, based on two-dimensional (2D)/three dimensional (3D) analytical solutions in finite and semi-infinite domains with zero-flux boundary conditions (BC). The RWPT algorithm induced from these analytical solutions is then applied to more complex geometries of grains and pores (that are difficult to solve analytically). Different configurations or structures at the microscale level (e.g., grains and pores in a porous media) will be chosen in order to obtain composite isotropic media with different porosities at the macroscale level. Then, by choosing elongated microstructures, anisotropic media emerge at the macroscopic level. These configurations will be obtained from microscale motifs of grains and pores with different porosities, repeated periodically and indefinitely (infinite domain grain–pore system). Effective macroscale properties (porosities, effective diffusion tensors, tortuosities) are then calculated for instantaneous injection (initial point source), using moments of particles positions. While standard homogenization models use steady state algorithms, the present RWPT allows the capture of the transient dynamics of the system. The method gives additional information on the evolution of the macro-properties through time, and thus ultimately, determines
the representative elementary volume (REV) of the system in addition to the macro-properties. In addition, RWPT is easily scalable with no limitations in the number of processors.

This paper is organized as follows. The next section, Section 2, presents the theory behind RWPT methods, and the corresponding concentration-based PDE, particularly for diffusion problems with zero-flux BC's. Section 3 presents a generalization of particle-based methods and algorithms for solving diffusion transport problems using RWPT in granular media with multiple zero-flux interfaces. In Section 4, the results obtained for different geometries of grains and pores are presented and discussed. An upscaling method using spatial moments (such as $\sigma_{xx}^2(t)$) is applied to the particle-based microscale simulation results in order to deduce macroscopic diffusion properties. The conclusive Section 5 recapitulates results and discusses extensions.

## 2. Theory: diffusion PDE, boundaries, and random particle positions

In this first part, we introduce the equations and the formalism in order to establish zero-flux conditions using particle methods.

### 2.1. PDE's for initial-value problems in infinite domains or with one zero-flux BC

#### 2.1.1. The Gaussian PDF of particle positions (infinite domain)

Let us define the initial-value problem (Cauchy's problem) in an infinite porous domain:

$$ egin{cases} 
\forall t > 0; \forall x \in \mathbb{R}; & \frac{\partial C}{\partial t}(x, t) = -V \frac{\partial C}{\partial x}(x, t) + D \frac{\partial^2 C}{\partial x^2}(x, t) \\ 
\forall t > 0; & \lim_{x \to \pm\infty} C(x, t) = 0 \\ 
\forall x \in \mathbb{R}; & C(x, 0) = \frac{M_0}{\theta} \delta(x - x_0), 
\end{cases} 
\quad (1) $$

where $\theta$ is the (constant) porosity, $M_0$ the initial mass of solute at $x = x_0$, and $C$ is the mass concentration of solute per volume of solvent (kg/m in 1D).

The solution of (1) is:

$$ \forall t > 0; \forall x \in \mathbb{R}; \quad C(x, t) = \frac{M_0}{\theta} \frac{1}{\sqrt{4\pi Dt}} \exp\left(-\frac{(x - (x_0 + Vt))^2}{4Dt}\right). \quad (2) $$

Equation (2) can be written as follows:

$$ \forall t > 0; \forall x \in \mathbb{R}; \quad C(x, t) = \frac{M_0}{\theta} G(x_0 + Vt, 2Dt; x), \quad (3) $$

where $G(x_0+Vt, 2Dt; x)$ is the Gaussian probability density function (PDF) with mean $x_0+Vt$ and variance $2Dt$. In particular, $G(0, 2t; x)$ is the fundamental solution of the diffusion equation [17].

$N(\mu, \sigma^2)$ denotes a Gaussian random variable (RV) with mean $\mu$ and variance $\sigma^2$, which means it has the PDF $x \to G(\mu, \sigma^2; x)$.

In 2D space, the initial-value PDE problem for pure diffusion ($V=0$) is:

$$ egin{cases} 
\forall t > 0; \forall (x, y) \in \mathbb{R}^2; & \frac{\partial C}{\partial t}(x, y, t) = D \left( \frac{\partial^2 C}{\partial x^2}(x, y, t) + \frac{\partial^2 C}{\partial y^2}(x, y, t) \right) \\ 
\forall t > 0; & \lim_{x,y \to \pm\infty} C(x, y, t) = 0 \\ 
\forall (x, y) \in \mathbb{R}^2; & C(x, y, 0) = \frac{M_0}{\theta} \delta(x - x_0, y - y_0). 
\end{cases} 
\quad (4) $$

The solution of (4) is:

$$ \forall t > 0; \forall x \in \mathbb{R}; \quad C(x, y, t) = \frac{M_0}{\theta} G((x_0, y_0), 2Dt; (x, y)), \quad (5) $$
where $G((x_0, y_0), 2Dt; (x, y))$ is the bi-variate Gaussian PDF for 2 i.i.d. variables (independent and identically distributed variables) with mean $(x_0, y_0)$ and isotropic variance $\sigma_X^2(t) = \sigma_Y^2(t) = 2Dt$. This bi-variate Gaussian solution of (4) will be used in subsequent subsections.

### 2.1.2. Semi-infinite domain with zero-flux BC

In a semi-infinite one-dimensional (1D) domain with a zero-flux BC at $x = x_{BC}$, the diffusion of an initial-value “source” located at $x = x_0$ is governed by the following system of equations:

$$ egin{cases} 
\forall t > 0; \forall x \in [-\infty; x_{BC}]; & \frac{\partial C}{\partial t}(x, t) = D \frac{\partial^2 C}{\partial x^2}(x, t) \\ 
\forall t > 0; & \lim_{x \to -\infty} C(x, t) = 0 \\ 
\forall t > 0; & -\theta D \frac{\partial C}{\partial x}(x_{BC}, t) = 0 \\ 
\forall x \in [-\infty; x_{BC}]; & C(x, 0) = \frac{M_0}{\theta} \delta(x - x_0). 
\end{cases} 
\quad (6) $$

The solution of (6) is (7) (similar to the one found in [18, p. 363–365]): $\forall t > 0; \forall x \in [-\infty; x_{BC}]$;

$$ C(x, t) = \frac{M_0}{\theta} (G(x_0, 2Dt; x) + G(2x_{BC} - x_0, 2Dt; x)). \quad (7) $$

Equation (7) could be written as

$$ C(x, t) = \frac{M_0}{\theta} (G(x_0, 2Dt; x) + L_{BC}(G(x_{BC}, 2Dt; x))), \quad (8) $$

where $L_{BC}$ is a linear application that transforms $G(x_0, 2Dt; x)$ into $G(2x_{BC} - x_0, 2Dt; x)$. Notice here that $L_{BC}(G(x_0, 2Dt; x))$ is the symmetric of $G(x_0, 2Dt; x)$ relative to $x = x_{BC}$ (see Figure 1).

In the 2D case considering the left half-space with a zero-flux BC on the vertical plane $x = x_{BC}$, Equation (6) becomes

$$ egin{cases} 
\forall t > 0; \forall x \in [-\infty; x_{BC}]; \forall y \in \mathbb{R}; & \frac{\partial C}{\partial t}(x, y, t) = D \left( \frac{\partial^2 C}{\partial x^2} + \frac{\partial^2 C}{\partial y^2} \right) (x, y, t) \\ 
\forall t > 0; & \lim_{x,y \to -\infty} C(x, y, t) = 0 \\ 
\forall t > 0; \forall x \in [-\infty; x_{BC}]; & \lim_{y \to +\infty} C(x, y, t) = 0 \\ 
\forall t > 0; \forall y \in \mathbb{R}; & -\theta D \frac{\partial C}{\partial x}(x_{BC}, y, t) = 0 \\ 
\forall x \in [-\infty; x_{BC}]; \forall y \in \mathbb{R}; & C(x, y, 0) = \frac{M_0}{\theta} \delta(x - x_0, y - y_0). 
\end{cases} 
\quad (9) $$

The solution of this equation (9) is the following combination of Gaussians: $\forall t > 0; \forall x \in [-\infty; x_{BC}]; \forall y \in \mathbb{R}$;

$$ C(x, y, t) = \frac{M_0}{\theta} (G((x_0, y_0), 2Dt; (x, y)) + G((2x_{BC} - x_0, y_0), 2Dt; (x, y))). \quad (10) $$

The half-domain solution (10) (with zero-flux BC) will be used in subsequent subsections.

### 2.2. PDE solutions with zero-flux BC's

The concentration-based analytical solutions shown in this section have been checked by direct substitution in their respective system of equations and their initial and BC's. They take the form of infinite series. Using the same concept shown in this section for 2D analytical solutions, these solutions could be extended to solve 3D problems. An analytical method, analogous to that described in Appendix A, will be applied later to RWPT for treating multiple interfaces in 2D (see RWPT algorithms in Section 3.2 further below).

The method used in Appendix A to obtain 2D concentration-based solutions $C(x, y, t)$ was applied to different initial-boundary-value problems with cartesian parallel interfaces serving as
<figure><img src="image_5.png" /><figcaption>**Figure 1.** Plot of different terms of (7) and (8) (analytical solutions). The red vertical line represents the interface: BC with zero flux. The green vertical line is the initial source position located at $x = 2$. The bold curves are the solutions inside the domain of interest, while the dashed lines are the same solutions extended outside the domain. The green bold curve is a Gaussian: it is the source solution in an infinite domain. The magenta Gaussian curve is the symmetric of the green one relative to the interface. The bold blue curve is the left-hand side of (7) and (8); it is the sum of the green and magenta curves.</figcaption></figure>

zero-flux planes. In these problems, each and every interface is infinite. Even when the domain is finite in one direction, the corresponding interface (orthogonal to that direction) is assumed to be infinite. This last condition restricts the application of this method to the geometric configurations of Figures 24–27. These configurations allow the solution to be the product of the 1D solution in one direction multiplied by the 1D solution in the other direction.

In the following Section 3.2, it is suggested that the independence between the two directions in these concentration-based solutions implies an independence between the $X(t)$ and $Y(t)$ positions of particles. This is one of the key points of our RWPT algorithms to be implemented for various grain-pore systems (see results in Section 4).

## 2.3. From particle positions to concentration

Let us consider now a particle-based method to solve initial-value problems. The concentration will be determined from the particles positions PDF (note, in the present case, total mass $M_0$ is constant, and concentration $C(x, t)$ is proportional to PDF):

$$C(x, t) = \int_{\mathbb{R}} C(X_t, t) \delta(X_t - x) dX_t \qquad (11)$$

$$C(x, t) = \int_{\mathbb{R}} \delta(X_t - x) dm_t, \qquad (12)$$
where $X_t$ and $dm_t$ represent, respectively, the position and mass of an infinitesimal concentration packet (to be discretized as a “particle”).

### 3. Methods and algorithms

Based on the previously established solution principles for the concentration-based PDE's (see Section 2.2), this section develops an RWPT algorithm that solves diffusion problems in more complex geometries (compared to the geometries of Section 2.2).

#### 3.1. RWPT methods with one zero-flux BC

RWPT methods are primarily inspired from initial-value problems (see Section 2.1). These methods consider particles to be independent. Therefore, if a particle does not encounter any interface, then it will be considered as if it were in an infinite domain, and follows the corresponding algorithm. On the other hand, the RWPT algorithm is altered for particles that do encounter interfaces (e.g., interface with zero-flux BC).

##### 3.1.1. RWPT in an infinite domain

To solve (1) using particle methods in 1D, for a time $t$, we need to generate particles with a PDF equivalent to (2). Thus, these particles must follow the RV $N(x_0 + Vt, 2Dt)$, hence

$$X_t = N(x_0 + Vt, 2Dt). \qquad (13)$$

Using simple Gaussian RV operations, we obtain

$$X_t = x_0 + Vt + \sqrt{2Dt}; N(0, 1), \qquad (14)$$

where $N(0, 1)$ is a Gaussian RV, with 0 mean and unit variance.

In conclusion, the PDF of RV $X_t$ is identical to the concentration solution in (2) divided by $M_0/\theta$.

In 2D, to solve a pure isotropic diffusion problem (4) using a particle method, $X_t$ governed by

$$\mathbf{X}_t = \begin{pmatrix} x_0 \\ y_0 \end{pmatrix} + \sqrt{2Dt}; \begin{pmatrix} N_X(0, 1) \\ N_Y(0, 1) \end{pmatrix}, \qquad (15)$$

where $N_X(0, 1)$ and $N_Y(0, 1)$ are two independent Gaussian RV's, with 0 mean and unit variance.

##### 3.1.2. RWPT in semi-infinite domain with zero-flux BC

Let us define $\Omega_{BC}$ and $\Omega_{BC}^*$ as: $\Omega_{BC} = ]-\infty; x_{BC}]$ and $\Omega_{BC}^* = [x_{BC}; +\infty[$, so $\Omega_{BC} \cup \Omega_{BC}^* = \mathbb{R}$. In Section 2.1.2, we have shown that the solution of the semi-infinite domain with zero-flux BC could be written as the sum of: the solution of the infinite domain $G$ and its symmetric $L_{BC}(G)$ relative to the position of the interface $x = x_{BC}$. Both are truncated on the semi-infinite domain $\Omega_{BC}$. On the other hand, $L_{BC}(G)$ truncated on $\Omega_{BC}$ is the symmetric of $G$ truncated on $\Omega_{BC}^*$. Therefore, in RWPT methods, $G$ truncated on $\Omega_{BC}$ represents the particles that did not reach the interface. The particles algorithm (i.e., infinite domain algorithm) is not altered; this is because $G$ is the solution of the infinite domain.<fnref n="1" /> Whereas, $L_{BC}(G)$ truncated on $\Omega_{BC}$ symmetry represents the particles that crossed the interface were totally reflected (i.e., their final position is the symmetric relative to the interface).

<footnote n="1">¹These considerations (properties) are due to the linearity of the corresponding PDE problem. A superposition principle is being applied; and there are analogies with Green's functions approach and with the method of images as well.</footnote>
Hence, the particles Random Walk algorithm becomes, in the presence of a zero-flux plane at $x = x_{BC}$:

$$ \left\{ \begin{array}{ll} \forall p & X_t^p = N^p(x_0, 2Dt) \\ \text{If } (X_t^p > x_{BC}) & X_t^p = 2x_{BC} - x_0 - \sqrt{2Dt}N^p(0, 1). \end{array} \right. \qquad (16) $$

Note here that $N^p$ is a Gaussian random variate, that is, a given realization of the Gaussian RV, for particle number $p$. Since its outcome is known, it can be compared to the interface position $x_{BC}$.

## 3.2. Multiple zero-flux interfaces for RWPT

In RWPT methods, after crossing one interface, particles may eventually cross multiple interfaces. First, because the spatial step depends on a Gaussian RV, which has unbounded values. Secondly, in multidimensional problems, particles located near a corner, for example, around two perpendicular interfaces, might cross the interfaces forming this corner [19–22]. Thus, the algorithm governing particles displacements must take into account such possibility. Therefore, this subsection proposes a generalization, analogous to the concept of Section 2.2, where the algorithm for a single interface with zero-flux boundary (16) was used to treat multiple interfaces.

### 3.2.1. RWPT zero-flux BC in 1D

Using an analogous method to the one used to solve concentration-based problem with two zero-flux BC's (see Appendix A.1), and combining it with the RWPT algorithm for one zero-flux BC described in Section 3.1.2, the 1D RWPT for two zero-flux BC's becomes Algorithm 1 below.

<caption>**Algorithm 1:** RWPT algorithm in a bounded domain with zero-flux BC's</caption>

$$ \left\{ \begin{array}{ll} \forall p & X^p(t^{n+1}) = X^p(t^n) + \sqrt{2D(t^{n+1} - t^n)} N^{n;p}(0, 1) \\ \textbf{While} & (X^p(t^{n+1}) \notin [x_{BC1}; x_{BC2}]) \\ \textbf{If} & (X^p(t^{n+1}) < x_{BC1}) \\ & X^p(t^{n+1}) = 2x_{BC1} - X^p(t^{n+1}) \\ \textbf{ElseIf} & (X^p(t^{n+1}) > x_{BC2}) \\ & X^p(t^{n+1}) = 2x_{BC2} - X^p(t^{n+1}) \\ \textbf{End} \\ \textbf{End} \end{array} \right. $$

Algorithm 1 calculates the final position of each particle using RWPT in a bounded 1D domain with zero-flux BC's. Each time a particle crosses an interface, its final position undergoes a total reflection.

The method used in Appendix A.1 to solve the concentration-based problem with two zero-flux BC's is similar to the use of the method of images through infinite superpositions. Each time the method of images was used, it corrected the new term to fit the boundary condition.

In RWPT, the method of images becomes a total particle reflection. If the new particle position is outside the domain $[x_{BC1}; x_{BC2}]$, it represents, in term of PDE's, the tail of the solution that it is outside of the domain and needs to be corrected again. That is why reflections are repeated until the particle's final position becomes inside the domain $[x_{BC1}; x_{BC2}]$. Note here that spatial steps lengths have a finite variance because they are Gaussian distributed. This algorithm (Algorithm 1) does not have restrictions on the time step. Therefore, the time step can be set equal to the final simulation time.
### 3.2.2. RWPT zero flux, interfaces in 2D

Section 2.2 presented 2D analytical solutions for various zero-flux configurations. In all these configurations, interfaces are infinite. This condition allows the concentration-based solutions to be the product of the 1D solution in one direction multiplied by the 1D solution in the other direction. In RWPT, if the PDE of the stochastic process could be written as the product of the PDE's of its components, then the position of X and Y are independent. In this case, there are no restrictions on the time step and the 2D RWPT algorithm can be expressed as Algorithm 1 applied to X and Y.

On the other hand, if interfaces are not infinite, then it can be demonstrated that the concentration-based solution cannot be written as a product of two functions with separate and independent variables. The proposed following algorithm (Algorithm 2) is a generalization of the infinite interface case taking into account the case of finite interfaces and complex geometries in a Cartesian mesh.

<caption>**Algorithm 2:** 2D RWPT for infinite or finite interfaces with zero-flux BC</caption>

(1) Consider particle (p) with mass $m_p$ and position $X_1^{(p)}$ and $X_2^{(p)}$.

(2) For $i \in \{1;2\}$; $X_i^{(p)}(t^{n+1}) = X_i^{(p)}(t^n) + \sqrt{2D(t^{n+1} - t^n)}N_i^{n;p}(0,1)$

(3) $\begin{cases} X_i^{(p)}(1) = X_i^{(p)}(t^{n+1}); & X_i^{(p)}(0) = X_i^{(p)}(t^n) \\ u = \operatorname{sign}(N_i^{n;p}(0,1)) & t^0 = 0 \end{cases}$

(4) **while** particle (k) crosses interfaces **do**

(a) Calculate the times $t_1$ and $t_2$ when (k) reach the nearest interface in the x and y directions

(b) $[I; t_I] = \min(t_1, t_2)$; $i$ is the direction that reached its interface first (I could be either direction 1 or 2) and $t_I$ is its corresponding time

(c) $\begin{cases} X_I^{(p)}(0) = X_{\text{Interface}} \\ X_{3-I}^{(p)}(0) = X_{3-I}^{(p)}(0) + u_{3-I}\sqrt{2D}(\sqrt{t_I} - \sqrt{t^0})N_I^{n;p}(0,1) \end{cases}$

(d) $\{u_I = -u_I \quad t^0 = t_I\}$

(e) $\{i \in \{1;2\}; \ X_i^{(p)}(1) = X_i^{(p)}(0) + u_i\sqrt{2D}(\sqrt{t^{n+1}} - t^n) - \sqrt{t_I}N_i^{n;p}(0,1)\}$

(5) **end**

(6) $X_i^{(p)}(t^{n+1}) = X_i^{(p)}(1)$

Algorithm 2 calculates the final position for particles using RWPT in a 2D domain with zero-flux interfaces. Each time a particle crosses an interface, its final position undergoes a 2D total reflection based on the analytical method of images. Particles may cross interfaces multiple times in the same time step. Every time they do, there their position is recalculated.

### 3.2.3. Discussion on time step size

Algorithm 1 requires particles to encounter interfaces that are only parallel to each other. Hence, the generalization to Algorithm 2. However, this algorithm (Algorithm 2) requires a time step small enough that most particles do not encounter multiple horizontal and vertical interfaces during a single time step. Therefore, the maximum spatial step of a particle should
<figure><img src="image_6.png" /><figcaption>**Figure 2.** Effective diffusivity versus time step. The bold blue curve represents the effective diffusivity plotted against the time step chosen for the simulation. The horizontal dashed red line correspond to the value of effective diffusivity using a time step $\Delta t = 10^{-3}$ (or roughly $\Delta t \le 10^{-2}$). The vertical dashed green line corresponds to the time step $\Delta t = 10^{-2}$ used in all subsequent simulations.</figcaption></figure>

not be greater than the characteristic length of the microgeometry. This limitation could be interpreted in terms of constraints as follows:

$$ \max(\Delta X_{\text{step}}) \le \lambda_{\text{PORE}} \qquad (17) $$

$$ \max(\sqrt{2D_{\text{PORE}}\Delta tZ}) \le \lambda_{\text{PORE}}, \qquad (18) $$

where Z is a normalized Gaussian RV $N(0, 1)$.

$$ \Delta t \le \frac{\lambda_{\text{PORE}}^2}{2D_{\text{PORE}}(\max(Z))^2}. \qquad (19) $$

For this problem of 2D grain/pore, the maximum absolute value of the Gaussian RV will not exceed 7 in practice (knowing that the probability of the absolute value of a normal Gaussian RV to exceed 7 is around $10^{-12}$, this probability is calculated from the cumulative distribution function of the normal distribution). Hence, we replace $(\max(Z))^2$ by 50, and (19) becomes

$$ \Delta t \le 10^{-2} \frac{\lambda_{\text{PORE}}^2}{D_{\text{PORE}}}. \qquad (20) $$

In addition, we should keep in mind that for the same simulation time, a small time step costs more in term of computational time compared to a larger time step. Hence, we choose the largest time step that verifies condition (20). In the following subsection, we take $\lambda_{\text{PORE}} = 1$ and $D_{\text{PORE}} = 1$, therefore the time step is equal to $\Delta t = 10^{-2}$. Furthermore, to validate $\Delta t = 10^{-2}$, we repeated the same simulations with different time steps. The result in terms of the calculated effective diffusion coefficient (explained later in Section 4.1.1) is shown in Figure 2 versus time step size. It is clear that the constraint $\Delta t \le 10^{-2}$ is adequate, and required.
<figure><img src="image_7.png" /><figcaption>**Figure 3.** 2D Micromodel motif for a 75% porosity. The red square represents the grain, the white region represents the pore.</figcaption></figure>

### 3.3. Grain-pore micromodel

#### 3.3.1. Grain/pore motif and pattern

This subsection describes the first geometric configuration that will be studied. Figure 3 represents the Grain/Pore motif that will be repeated, periodically and indefinitely, to create the infinite domain geometry. The red square corresponds to the grain while the pore is in white. The area of the pore divided by the area of the whole motif (grain/pore), is the porosity of the motif (in this case $\theta = 75\%$). Since the infinite domain geometry is a periodic repetition of the previous motif (see Figure 4), then the effective porosity of the infinite domain is the same: $\theta_{eff} = 75\%$.

#### 3.3.2. Representative elementary volume (REV)

For this kind of periodic structure, we could naively imagine that the REV would be composed of one grain and three pores as shown in Figure 3. However, for this particular problem of the diffusion from an initial source in an infinite domain, we will later define the REV using a convergence criterion for macroscopic properties (such as effective diffusivity).

The first macroscopic property to be studied here is the porosity ($\theta$). To define the REV using the macroscopic porosity, we define first a square surface (in 3D, we would have used a volume) in which the effective porosity is calculated. In most cases, the surface is taken as a disc (sphere in 3D), but here, in the interest of simplicity, we take a square. This area will be related to a “radius” (distance from the center (the initial source position) to the area’s edge). This “radius” is equal to half the edge of the square, in the case of a disc, this distance is the radius of the disc. Later we will use this “radius” to define the limits of our REV.

Figure 5 plots the porosity inside a square against the “radius” of this square. In Figure 4, the blue dot represents the center of the square, and the three magenta squares are examples of areas over which the porosity is calculated. The three magenta vertical lines in Figure 5 correspond to the magenta squares outlined in Figure 4.

**Simulation time.** The simulation time should be large enough for the particles to explore the domain consisting of one or several representative elementary volumes (REVs). Thus, the reference size specific to the simulation time ($L_{REF}$) could be a few times the typical REV length scale ($L_{REV}$).
<figure><img src="image_8.png" /><figcaption>**Figure 4.** 2D Micromodel comprising 20 × 20 grains (red squares). The blue dot in the center is used as the initial source. The three square frames in magenta represent upscaling domains of various sizes.</figcaption></figure>

The simulation time can be estimated by equating the reference size ($L_{REF}$) to the standard deviation displacement (i.e., the so-called root-mean-square dispersion length). This yields the estimated simulation time for the particle cloud to explore a few times the typical REV length scale of the grain/pore system:

$$T_{\text{Simu}} \sim \frac{L_{\text{REF}}^2}{2D_{\text{eff}}}. \qquad (21)$$

# 4. Results and discussion

This section analyzes the results of simulations of diffusion of an initial source inside the pores of a composite media pore/grain. Different structures are studied, either isotropic or anisotropic and with different porosities. For each case, an upscaling method, using the time-dependent moments of particles positions, is applied to determine the macro-properties of the composite media (essentially, effective porosity tortuosity, and diffusion coefficient).

## 4.1. Isotropic grains ($\lambda_{\text{Grain},X} = \lambda_{\text{Grain},Y}$)

### 4.1.1. Porosity $\theta = 75\%$

The first grain/pore configuration, where the diffusion (Algorithm 2) will be tested, is the one described in Figure 4. An initial source diffuses inside the center of the pore at the position (0.5; 0.5). The result of the simulation is shown in Figure 6. The blue dots represents the particles
<figure><img src="image_9.png" /><figcaption>**Figure 5.** Porosity profile versus “radius” of the upscaling domain expressed in pore size units in bold blue curve. The horizontal bold red line is the effective porosity of the infinite domain 75%. The horizontal dashed red lines are ±5% of the bold red line. The vertical magenta bold lines correspond to the three magenta squares outlined in Figure 4.</figcaption></figure>

<figure><img src="image_10.png" /><figcaption>**Figure 6.** 2D plot of particles positions (blue dots): result of diffusion simulation at time $t_F > T_{\text{Simu}}$ in the geometry of Figure 4.</figcaption></figure>
<figure><img src="image_11.png" /><figcaption>**Figure 7.** Zoom at the center of Figure 6: 2D plot of particles positions (blue dots) results of diffusion simulation at time $t_F > T_{Simu}$, in the geometry of Figure 4.</figcaption></figure>

after time $t_F = 20$. This time is chosen so that enough particles have explored the periodic infinite domain, which means that the particles spread over more than one REV.

Figure 7 is a zoom of Figure 6 centered on the position of the initial source (0.5; 0.5). It shows clearly that particles are found anywhere but inside the squared grains. Algorithm 2 makes particles bounce off the grains.

After post-processing, as described in Section 2.3, Figure 8 is obtained. In this figure, concentrations $C(x, y)$ were computed as 2D histograms based on the smallest histogram bins corresponding pore size $\lambda_{PORE}$. They represent the 2D concentrations of particles. As one might predict, the highest concentrations are found near the center where the initial source was (0.5; 0.5). Then, the concentration drops gradually while moving further away from the center, as expected.

**Effective diffusivity.** There are different methods to calculate the effective diffusivity. The most common one in the literature uses the concentration gradient in a Dirichlet bounded domain [23]. This study uses the second order spatial moments of particles positions $M_2(t)$:

$$M_2(t) = \begin{pmatrix} \sigma_X^2 & \sigma_{XY} \\ \sigma_{YX} & \sigma_Y^2 \end{pmatrix}. \qquad (22)$$

If $M_2(t)$ is diagonal then $M_2(t) = \sigma^2(t)I$, with $I$ is the identity matrix. Then, the effective diffusivity is determined by one of the two following equations:

$$D_{\text{eff}} = \frac{\sigma^2}{2t} \qquad (23)$$

$$D_{\text{eff}} = \frac{1}{2} \frac{d\sigma^2}{dt}. \qquad (24)$$

Figures 10 and 11 plot the effective diffusivity as defined in (23). The convergence of the effective diffusivity is defined as: it is convergent when its value stays between ±1% of its final mean value. This happens at around $t = 1$ close to the value predicted by (21).
<figure><img src="image_12.png" /><figcaption>**Figure 8.** 2D plot of concentrations of the simulated particle cloud shown in Figure 6.</figcaption></figure>

Initially, the effective diffusivity is equal to the local diffusion coefficient inside the pore. $D_{\text{eff}}$ stays equal to $D_0$ until particles start encountering grains interfaces. Then, as times goes by, particles bounce off the grains, which limits their diffusion, and $D_{\text{eff}}$ starts slowly to decline, until it converges toward its final value (because of the periodicity of the geometry), which is smaller than $D_0$ due to tortuosity effect.

The REV can be defined from the convergence in time of macro-properties. From Figure 11, we determine that effective diffusivity converges at time equal to 1. Hence, the effective diffusivity REV would have a “radius” of about $\sqrt{2D_{\text{eff}}t}\max(Z)$. Therefore, the effective diffusivity REV “radius” is around 8.5. On the other hand, the geometrical porosity REV “radius” is around 7 (see Figure 5). The two REV's sizes are relatively close, which could be explained by the dependency between the porosity and the effective diffusivity (and their relative convergence). The small difference between the two REV's sizes could be explained by the approximation shape of the REV's, square for the porosity REV and circle for the effective diffusivity REV.

Figure 9 shows the trajectories of 3 particles. The particles circumvent the grains (red squares), and that is what creates the tortuosity of the system.

As a comparison, we will now vary the porosity of the medium in order to investigate the effects of porosity on the effective diffusivity.

### 4.1.2. *Experiments with different porosities*

Table 1 shows the results of different simulations with different Grain/Pore configurations. In all the configurations tested in this table, the grains are taken as squares, characterized by their size $\lambda_{\text{Grain}}$. Usually, we refer the tortuosity to be the ratio of the molecular diffusion to the effective diffusion value. In the literature, we can find several definitions and empirical laws
<figure><img src="image_13.png" /><figcaption>**Figure 9.** Plot of the trajectories (green, brown, and magenta) of 3 particles inside the grain/pore system of Figure 4: the red squares are the grains and the white regions are the pores. The bold blue dot is the initial position $x = (0.5; 0.5)$ of the three particles. The green, brown, and magenta stars are the final positions of the particles.</figcaption></figure>

<figure><img src="image_14.png" /><figcaption>**Figure 10.** Effective diffusivity versus time. The blue and green bold curves are the variance of the positions X and Y respectively divided by $(2 * t)$ with t is the time. The red bold line is the mean value of the variance 0.753. The red dashed lines are ±1% of the bold red line.</figcaption></figure>
<figure><img src="image_15.png" /><figcaption>**Figure 11.** Log-log of Figure 2 effective diffusivity versus time. The blue and green bold curves are the variance of the positions *X* and *Y* respectively divided by (2 * t) with *t* is the time. The red bold line is the mean value of the variance 0.753. The red dashed lines are ±1% of the bold red line.</figcaption></figure>

<table><caption>**Table 1.** Porosity versus effective diffusivity versus tortuosity for different configurations</caption><thead><tr><th>λ<sub>Grain</sub>/λ<sub>Pore</sub></th><th>1</th><th>4</th><th>9</th><th>19</th><th>49</th><th>99</th></tr></thead><tbody><tr><th>θ</th><td>75%</td><td>36%</td><td>19%</td><td>9.75%</td><td>3.96%</td><td>1.99%</td></tr><tr><th>D<sub>eff</sub></th><td>0.753</td><td>0.586</td><td>0.542</td><td>0.521</td><td>0.508</td><td>0.504</td></tr><tr><th>τ<sup>2</sup></th><td>1.33</td><td>1.71</td><td>1.84</td><td>1.92</td><td>1.97</td><td>1.98</td></tr><tr><th>θD<sub>eff</sub>, from this work</th><td>0.565</td><td>0.211</td><td>0.103</td><td>0.051</td><td>0.02</td><td>0.01</td></tr><tr><th>θD<sub>eff</sub> from [26]</th><td>0.583</td><td>0.213</td><td>—</td><td>—</td><td>—</td><td>—</td></tr></tbody></table>

used to determine tortuosity values. For instance, [24,25] have developed correlations in order to determine the tortuosity values, they are directly related to the gas saturation and the porosity.

Here the porosity is calculated using the formula (25):

$$ \theta = \frac{1 + \lambda_{\text{Grain},X} + \lambda_{\text{Grain},Y}}{1 + \lambda_{\text{Grain},X} + \lambda_{\text{Grain},Y} + \lambda_{\text{Grain},X} + \lambda_{\text{Grain},Y}}. \quad (25) $$

With isotropic grains we have $λ_{Grain,X} = λ_{Grain,Y} = λ_{Grain}$, and $λ_{Pore} = 1$.

The tortuosity is calculated according to the definition of [23]:

$$ \tau^2 = \frac{D_m}{D_{\text{eff}}} \quad (26) $$

with $D_m$ is the diffusion coefficient inside the pore (also denoted here $D_0$).

For the two first columns, we are in an excellent agreement with the results reported by [26]. Unfortunately, the results for the other porosities are not available in order to completely validate this set of numerical experiments.

Thus, the results of this work are compared with the following models (plotted in Figure 12):

* Maxwell [27] analyzed analytically a dilute suspension of spheres and obtained:

$$ \theta D_{\text{eff}} = \frac{\theta}{1 + (1 - \theta)/2}. \quad (27) $$

* Weissberg [28] used a variational method and proposed the following:

$$ \theta D_{\text{eff}} = \frac{\theta}{1 - \ln(\theta)/2}. \quad (28) $$
<figure><img src="image_16.png" /><figcaption>**Figure 12.** Plot of $\theta D_{eff}$ versus $\theta$. The black, magenta, green, and red curves correspond, respectively, to $\theta D_{eff}$ of Maxwel, Weissberg, $\theta^{1.4}$ and $\theta^2$ models. The green ×, brown + correspond, respectively, to Ryan's and Kim's work. The bold blue circles correspond to the results of this work in Table 1.</figcaption></figure>

* Wakao and Smith [29] using a micropore macropore model proposed:

$$ \theta D_{eff} = \theta^2. \qquad (29) $$

* Kim et al. [30] proposed the empirical formulation:

$$ \theta D_{eff} = \theta^{1.4}. \qquad (30) $$

Figure 12 shows that the results of this study (the blue circles) are relatively close to results found in the literature. The small differences might be due to the difference of grain shapes (e.g., Maxwell [27] uses spherical grains for his study), or due to the grains arrangements (which are different between studies [28–30]). However, for this particular isotropic periodic case, the analytical expressions give very good approximations of the effective diffusivity.

Figure 13 shows clearly that the effective tortuosity of the Grain/Pore geometry tends toward the geometric tortuosity:

$$ \tau = \frac{\lambda_{Grain,X} + \lambda_{Grain,Y}}{\sqrt{\lambda_{Grain,X}^2 + \lambda_{Grain,Y}^2}} \qquad (31) $$

as found in the literature [31]. When $\lambda_{Grain,X} = \lambda_{Grain,Y}$, then $\tau = \sqrt{2}$.

Figure 14 plots the tortuosity versus the porosity. This figure shows that for $\theta \in [0;0.36]$, the tortuosity is linear with the porosity with a determination coefficient $R^2 = 99.99\%$ and the corresponding empirical formula is

$$ \tau \approx \sqrt{2} - 0.3\theta = \tau_{geo} - 0.3\theta. \qquad (32) $$

## 4.2. Anisotropic grains ($\lambda_x > \lambda_y$)

This subsection studies micro- and macroscopic diffusion in anisotropic geometries.
<figure><img src="image_17.png" /><figcaption>**Figure 13.** Plot of Table 1. Red, blue, and green bold curves correspond, respectively, to the squared tortuosity, the effective diffusivity, and the porosity plotted versus the anisotropy ratio. The dashed red, blue, and green lines correspond to the squared geometric tortuosity ($\tau_{geo}^2 = 2$), effective diffusivity ($D_{eff;geo} = 0.5$), and porosity ($\theta_{geo} = 0$).</figcaption></figure>

<figure><img src="image_18.png" /><figcaption>**Figure 14.** Tortuosity versus porosity. The bold blue curve with diamonds is the tortuosity from $\theta \in [0; 0.36]$. The green line is a linear ($\tau = 1.41 - 0.3\theta$) fitting of the blue curve with an $R^2 = 99.99\%$. The red curve with circles is the tortuosity from $\theta \in [0; 0.75]$.</figcaption></figure>

### 4.2.1. Representative elementary volume (REV) for anisotropic grains

Similar to the isotropic structure, the anisotropic REV is larger than one grain and three pores. Particles need to explore this structure multiple times to extract macroscopic properties. The REV is defined from the convergence of macroscopic properties. For elongated grains, the REV is elongated in the same direction as the grains. On the other hand, since particles X and Y displacements are dependent, X and Y macroscopic properties reach their convergence at the same simulation time.
<figure><img src="image_19.png" /><figcaption>**Figure 15.** 2D concentrations of particles in anisotropic geometry with $\theta = 75\%$. The ratio between $\lambda_x$ and $\lambda_y$ is 8.</figcaption></figure>

### 4.2.2. Anisotropic geometry with porosity $\theta = 75\%$

Figure 15 shows 2D concentration histogram of the result of an anisotropic test case with elongated grain $\lambda_x/\lambda_y = 8$, and with the same porosity as previously used $\theta = 75\%$. According to the previous analysis, we can extract from our results the values of $D_{eff;X} = 0.64$ and $D_{eff;Y} = 0.5$. Then it could be compared with the results available in [30]. In a similar case, with $\lambda_x/\lambda_y = 7.5$, they found a value of $D_{eff;X} = 0.717$ and $D_{eff;Y} = 0.175$. However, it might be noticed that the microstructure is not exactly the same as the solid particles are staggered and not in line, in addition, the shape of the pores in [30] are elongated $\lambda_{XPore}/\lambda_{YPore} = 0.77$ while the shape in the example of this study is taken as $\lambda_{XPore}/\lambda_{YPore} = 8$. So, it could generate different values of tortuosity.

### 4.2.3. Anisotropic geometry with porosity $\theta = 36\%$

**Anisotropy ratio** $\lambda_x/\lambda_y = 2$. Figures 16–18 shows the results of the diffusion of an initial source inside the center of a pore. The domain is a periodic repetition of a motif grain/pore with rectangular grains. The grain's length is two times longer than its width, and six times longer than the width of the pore. The periodic repetition of this motif creates an anisotropic infinite domain with a porosity of $\theta = 36\%$.

Figure 19 and its zoom Figure 20 plots the mean normalized concentration on each direction and divided by the local porosity to obtain the concentration inside the pores. These 1D effective
<figure><img src="image_20.png" /><figcaption>**Figure 16.** 3D Histogram of particles in anisotropic 2D geometry with $\theta = 36\%$. The ratio between $\lambda_x$ and $\lambda_y$ is 2.</figcaption></figure>

<figure><img src="image_21.png" /><figcaption>**Figure 17.** 3D Contour of particles positions in anisotropic 2D geometry with $\theta = 36\%$. The ratio between $\lambda_x$ and $\lambda_y$ is 2.</figcaption></figure>

concentration profiles were obtained by directional averaging in x and y:

$$ \overline{C^y}(x) = \frac{1}{Ly} \int C(x, y) dy. \qquad (33) $$
<figure><img src="image_22.png" /><figcaption>**Figure 18.** 2D concentrations of particles in anisotropic 2D geometry with $\theta = 36\%$. The ratio between $\lambda_x$ and $\lambda_y$ is 2.</figcaption></figure>

<table><caption>**Table 2.** Macroscopic properties (effective diffusion and tortuosity on the x and y directions) for different configurations (different ratios of lengths of elongated grains) with a fixed porosity $\theta = 36\%$.</caption><thead><tr><th>λ<sub>x</sub>/λ<sub>y</sub></th><th>4/4</th><th>6/3</th><th>24/2</th></tr></thead><tbody><tr><th>θ</th><td>36%</td><td>36%</td><td>36%</td></tr><tr><th>D<sub>eff;X</sub></th><td>0.586</td><td>0.726</td><td>0.930</td></tr><tr><th>D<sub>eff;Y</sub></th><td>0.586</td><td>0.429</td><td>0.126</td></tr><tr><th>τ<sup>2</sup><sub>x</sub></th><td>1.71</td><td>1.38</td><td>1.08</td></tr><tr><th>τ<sup>2</sup><sub>y</sub></th><td>1.71</td><td>2.33</td><td>7.94</td></tr></tbody></table>

It is clear that the diffusion inside the pores at the microscopic level (inside the pores) behaves as a pure diffusion in the macroscopic level with the macro-properties calculated using moments particles positions (see Figure 21).

Table 2 displays the results of different simulations in different configurations of grains/pores with the same porosity $\theta = 36\%$ and different anisotropy ratios $a = \lambda_x/\lambda_y$. For any given porosity $\theta_0$, and any given anisotropy ratio “a”, Equation (34) gives $(\lambda_x, \lambda_y)$ such that $\lambda_x/\lambda_y = a$, and such that the desired macroscale porosity $\theta_0$ is satisfied:

$$ \left\{ \begin{array}{l} \lambda_y = \frac{1 - \theta_0}{\theta_0} + \frac{a}{2a} \left( 1 + \sqrt{1 + \frac{4a}{(1+a)^2} \frac{\theta_0}{1 - \theta_0}} \right) \\ \lambda_x = a\lambda_y. \end{array} \right. \qquad (34) $$

### 4.2.4. Anisotropic geometry with porosity $\theta = 19\%$ and $\theta = 9.75\%$

The results shown in Tables 2–4 suggest the following trends:

* The longer a grain is in one direction, the higher its effective diffusion is in this same direction.
<figure><img src="image_23.png" /><figcaption>**Figure 19.** Effective concentration profiles C(x) and C(y) for anisotropic microgeometry: the bold blue and red curves correspond, respectively, to the mean (on direction x and y) concentration plotted against direction x and y. The bold green and magenta curves correspond, respectively, to the analytical Gaussian functions with diffusive coefficient $D_{eff;X} = 0.726$ and $D_{eff;Y} = 0.429$.</figcaption></figure>

* When the anisotropic ratio $\lambda_x/\lambda_y$ tends toward infinity, the effective diffusion tends toward 1 in the direction of the elongated grains (X), while in the opposite direction (Y) it tends toward 0.

* For a constant porosity, we have observed that the sum of the effective diffusivity of both directions stays constant regardless of the anisotropic ratio (see Figure 22).

The plot of effective diffusivity versus anisotropy ratio in Figure 22 shows that: for a constant anisotropy ratio, effective diffusivities of both directions increase with the porosity. In addition, for the configuration studied in this paper (rectangular periodic grains), the effective diffusivity components along the X and Y directions are related to the isotropic effective diffusivity through the following equation

$$D_{\text{eff;iso}} = \frac{D_{\text{eff;X}} + D_{\text{eff;Y}}}{2} \qquad (35)$$

Equation (35) holds for a given (fixed) porosity. $D_{eff,X}$ and $D_{eff,Y}$ are effective diffusivity components along the X and Y directions for the same anisotropic configuration, while $D_{eff,iso}$ is the effective diffusivity of the isotropic configuration ($\lambda_x/\lambda_y = 1$).

Figure 23 shows a fit of the tortuosity obtained numerically with the following analytical function:

$$\tau_{\text{fit;Y}} = (\tau_{\text{iso}} - \tau_{\infty}) \left(\frac{\lambda_x}{\lambda_y}\right)^{(3/4)} + \tau_{\infty} \qquad (36)$$

with $\tau_{iso}$ is the tortuosity in an isotropic configuration and $\tau_{\infty} = 1$ is the X tortuosity when $\lambda_x/\lambda_y \to \infty$.
<figure><img src="image_24.png" /><figcaption>**Figure 20.** Effective concentration profiles C(x) and C(y): zoom on the center of Figure 19.</figcaption></figure>

<figure><img src="image_25.png" /><figcaption>**Figure 21.** 2D concentrations of particles in anisotropic geometry with $\theta = 36\%$. The ratio between $\lambda_x$ and $\lambda_y$ is 12.</figcaption></figure>
<figure><img src="image_26.png" /><figcaption>**Figure 22.** Effective diffusivity versus anisotropy ratio $\lambda_x/\lambda_y$: the bold green curves with stars are for porosity $\theta = 36\%$. The bold red curves with squares are for porosity $\theta = 19\%$. The bold blue curves with circles are for porosity $\theta = 9.75\%$. For each color, the higher curve is $D_{\text{eff},X}$ and the lower curve is $D_{\text{eff},Y}$.</figcaption></figure>

<figure><img src="image_27.png" /><figcaption>**Figure 23.** Tortuosity versus anisotropy ratio $\lambda_x/\lambda_y$ for a porosity $\theta = 9.75\%$. The blue bold curve with circles is $\tau_x$. The bold red circles are $\tau_y$. The bold green curve is the analytical function: $\tau_{\text{fit};Y} = (\tau_{\text{iso}} - \tau_{\infty})(\lambda_x/\lambda_y)^{(3/4)} + \tau_{\infty}$.</figcaption></figure>

This fit (36) is in agreement with theoretical expectations for extreme anisotropy ratio $\lambda_x/\lambda_y \to \infty$ (diffusion occurs only on the X direction), and for the isotropic case $\lambda_x = \lambda_y$, where $\tau_{\text{fit};Y} = \tau_{\text{iso}}$.
<table><caption>**Table 3.** Macroscopic properties (effective diffusion and tortuosity on the *x* and *y* directions) for different configurations (different ratios of lengths of elongated grains) with a fixed porosity $\theta = 19\%$</caption><thead><tr><th>λ<sub>x</sub>/λ<sub>y</sub></th><th>9/9</th><th>10/8</th><th>12/7</th><th>17/6</th><th>34/5</th></tr></thead><tbody><tr><th>θ</th><td>19%</td><td>19%</td><td>19%</td><td>19%</td><td>19%</td></tr><tr><th>D<sub>eff;X</sub></th><td>0.542</td><td>0.593</td><td>0.654</td><td>0.763</td><td>0.881</td></tr><tr><th>D<sub>eff;Y</sub></th><td>0.542</td><td>0.500</td><td>0.412</td><td>0.309</td><td>0.159</td></tr><tr><th>τ<sup>2</sup><sub>X</sub></th><td>1.84</td><td>1.686</td><td>1.53</td><td>1.31</td><td>1.13</td></tr><tr><th>τ<sup>2</sup><sub>Y</sub></th><td>1.84</td><td>2.00</td><td>2.40</td><td>3.23</td><td>6.27</td></tr></tbody></table>

<table><caption>**Table 4.** Macroscopic properties (effective diffusion and tortuosity on the *x* and *y* directions) for different configurations (different ratios of lengths of elongated grains) with a fixed porosity $\theta = 9.75\%$</caption><thead><tr><th>λ<sub>x</sub>/λ<sub>y</sub></th><th>19/19</th><th>20/18</th><th>22/17</th><th>23/16</th><th>26/15</th><th>29/14</th><th>35/13</th><th>44/12</th><th>64/11</th><th>137/10</th></tr></thead><tbody><tr><th>θ</th><td>9.75%</td><td>9.75%</td><td>9.75%</td><td>9.75%</td><td>9.75%</td><td>9.75%</td><td>9.75%</td><td>9.75%</td><td>9.75%</td><td>9.75%</td></tr><tr><th>D<sub>eff;X</sub></th><td>0.521</td><td>0.534</td><td>0.583</td><td>0.605</td><td>0.628</td><td>0.673</td><td>0.715</td><td>0.782</td><td>0.852</td><td>0.948</td></tr><tr><th>D<sub>eff;Y</sub></th><td>0.521</td><td>0.512</td><td>0.467</td><td>0.431</td><td>0.396</td><td>0.360</td><td>0.293</td><td>0.235</td><td>0.164</td><td>0.078</td></tr></tbody></table>

# 5. Conclusion and perspectives

## 5.1. Conclusion

This study proposed an RWPT algorithm based on 2D/3D analytical solutions (that we have verified and proved the convergence) for zero-flux conditions. The new algorithm handles multiple reflections in the same time step, in addition to making the time step larger with less restrictions. With these larger time steps, the RWPT algorithm becomes computationally more efficient compared to previous algorithms. This efficiency is more prominent for anisotropic configurations. The algorithm was applied for diffusion of an initial source in an infinite domain (i.e., domain large enough so that the particles do not reach the domain edges during the simulation) with different configurations of grains and pores (different porosities in isotropic systems, and constant porosity with different anisotropic configurations). The aim is a multi-scale study from the microscale of grains and pores to the macroscale level. It leads to macro-properties (porosity, effective diffusivity, tortuosity in different directions). The upscaling method uses time-dependent moments of particles positions. While previous studies used steady state simulations with periodic boundary conditions. In this study, we used dynamic transient simulations in an infinite domain, which ultimately allowed us to define the REV from the convergence in time of macro- properties. The time step sensitivity study validated the time step choice for simulations. Results show the scale dependence of effective diffusion dynamically through both space and time scales. And the macro-properties converge toward the geometric theoretical values when the Grain/Pore geometry's porosity tends toward zero. Finally, results lead to new correlations that can be used to better understand the upscaling of media properties.

## 5.2. Perspectives

The strong analogy between the stochastic and the deterministic approaches allowed us to propose a new 2D RWPT algorithm. Similarly, using 3D deterministic approaches, the model
can be extended to simulation in 3D instead of 2D. The new less restrictive criterion on the time step can also help reduce the calculation time. Simulation results obtained not only complete previous upscaling results, but also extrapolate them into new empirical correlations that can help expand our knowledge of the upscaling. Moreover, a velocity component could be added to the RWPT algorithm, either by combining advective and diffusive transport in a single step, or by separating them with one diffusive sub-step followed by one advective sub-step. In addition, the grains can be replaced by porous aggregates with their own local porosity and diffusion coefficient (e.g., solute diffusion in a dual porous medium, or possibly, pressure diffusion in a dual porous medium). And finally, the deterministic microgeometries of pore-grain or pore-aggregate structures could be randomized in order to study their macroscale properties in terms of their random structure using the same random walk particle-based approach.

## Acronyms used in this text

* BC Boundary conditions

* PDE Partial differential equations

* PDF Probability density function

* 1D One-dimensional

* 2D Two-dimensional

* RWPT Random walk particle tracking

* REV Representative elementary volume

# Appendix A. Detailed PDE solutions under various geometries with zero-flux BC's

## A.1. 1D zero flux, finite domain

This subsection describes the method to solve 1D zero-flux problems in finite domain. This method will be used further below to propose the RWPT algorithm for multiple zero-flux BC.

As explained in Sections 2.1.2 and 3.1.2, the solution of a semi-infinite domain with a zero-flux boundary condition can be expressed as: the sum of the solution in an infinite domain G and of its symmetric relative to the position of the boundary $L_{BC1}(G)$. This superposition approach takes advantage of the linearity of the diffusion PDE, and it is analogous to the method of images. Thus, $L_{BC1}(G)$ added to G could be considered as a correction so that it fits the first zero-flux condition (BC1), noting that the second condition BC2 is already taken into account in the analytical solution "G" of (3).

$$ egin{cases} 
\forall t > 0; \forall x \in [x_{BC1}; x_{BC2}]; & \frac{\partial C}{\partial t}(x, t) = D \frac{\partial^2 C}{\partial x^2}(x, t) \\ 
\forall t > 0; & -\theta D \frac{\partial C}{\partial x}(x_{BC1}, t) = 0 \\ 
\forall t > 0; & -\theta D \frac{\partial C}{\partial x}(x_{BC2}, t) = 0 \\ 
\forall x \in [x_{BC1}; x_{BC2}]; & C(x, 0) = \frac{M_0}{\theta} \delta(x - x_0). 
\end{cases} 
\quad (37) $$

Therefore, in this paper, to solve the problem of (37) with two zero-flux BC's, we propose to apply the same concept for both BC's, that is, adding also $L_{BC2}(G)$ to the solution G so that it fits zero-flux BC number 2 too. However, although G has now been corrected to fit both BC's, $L_{BC1}(G)$ and $L_{BC2}(G)$ do not fit zero-flux BC2 and BC1, respectively. Hence, once again, using the same algorithm $L_{BC1}(G)$ is corrected by adding $L_{BC2}(L_{BC1}(G))$, so that $L_{BC1}(G) + L_{BC2}(L_{BC1}(G))$ does fit zero-flux BC2. Similarly, $L_{BC2}(G)$ is corrected by adding $L_{BC1}(L_{BC2}(G))$, so that $L_{BC2}(G) +
$L_{BC1}(L_{BC2}(G))$ does fit zero-flux BC1. Equation (38) shows how G (the solution of the infinite domain) is corrected by adding the terms to fit BC1 and BC2, respectively. Equation (38) shows that for each correction, we add a new layer of corrections ad infinitum:

$$C = G + L_{BC1}(G) + L_{BC2}(G) + L_{BC2}(L_{BC1}(G)) + L_{BC1}(L_{BC2}(G)) + \dots \quad (38)$$

Equation (38) becomes an infinite series (39) after applying the concept of correction (or super-position) indefinitely.

$$C = [i = 0] + \sum (\mathrm{id} + L_{BC1})(L_{BC2}L_{BC1})^i (\mathrm{id} + L_{BC2})(G). \quad (39)$$

Equation (39) can also be written as the following infinite series:

$$\begin{aligned} \forall t > 0; \quad & \forall x_{BC1} \le x \le x_{BC2}; \quad \frac{\theta}{M_0} C_{BC12}^x(x, t) \\ &= G(x_0, 2Dt, x) + [i = 1] + \sum (G(-2il_{x1-2} + x_0, 2Dt, x) \\ & \quad + G(2il_{x1-2} + x_0, 2Dt, x) + G(-2il_{x1-2} + x_{BC1} - x_0, 2Dt, x) \\ & \quad + G(2il_{x1-2} + x_{BC2} - x_0, 2Dt, x)), \end{aligned} \quad (40)$$

where $l_{x1-2} = x_{BC1} - x_{BC2}$.

## A.2. 2D zero flux, finite in x and infinite in y

This subsection describes the solution for a 2D zero-flux problem, finite in x and infinite in y directions. This configuration (see Figure 24) correspond to a particle positioned between two infinite planes, it will be used to propose the 2D RWPT algorithm for multiple zero-flux boundary conditions.

In 2D, the solution is the product obtained by multiplying the 1D solution in x direction and in y direction. Thus, the solution in the x and y directions are independent. In Figure 24, the 2D zero-flux finite domain in x direction and infinite in y direction, there is no BC in the y direction. Therefore, the solution is a Gaussian in the y direction multiplied by the solution in the x direction (40) obtained in Appendix A.1. This yields

$$\forall t \ge 0; \quad \forall x_{1:2} \le x \le x_{2:3}; \quad \forall y \in \mathbb{R};$$

$$\frac{\theta}{M_0} C(x, y, t) = G(y_0, 2Dt, y) \times \frac{\theta^x}{M_0} C_{BC12}^x(x, t). \quad (41)$$

## A.3. 2D zero flux, semi-infinite domain in x and y directions

The configuration studied in this subsection (Figure 25) correspond to a case where a particle is located near the intersection of two interfaces.

In this problem, there is one zero-flux BC for the x and y directions. Hence, the principle of correction (explained in Appendix A.1) must be applied so that the solution fits the zero-flux BC in x and y directions (instead of BC1 and BC2 of the same direction as in Appendix A.1). Since x and y are independent, then the solution of this problem (see Figure 25) is the solution of the 1D semi-infinite problem in x multiplied by the solution of the 1D semi-infinite problem in y. Hence, in this case the solution is not an infinite series but only the sum of four Gaussian functions.

$$C = (\mathrm{id} + L_{BC}^x)(\mathrm{id} + L_{BC}^y)(G) \quad (42)$$

$$\frac{\theta}{M_0} C(x, y, t) = G((x_0, y_0), 2Dt, (x, y)) + G((2x_{12} - x_0, y_0), 2Dt, (x, y)) \\ + G((x_0, 2y_{23} - y_0), 2Dt, (x, y)) + G((2x_{12} - x_0, 2y_{23} - y_0), 2Dt, (x, y)). \quad (43)$$
<figure><img src="image_28.png" /><figcaption>Figure 24. 2D parallel zero-flux interfaces.</figcaption></figure>

<figure><img src="image_29.png" /><figcaption>Figure 25. 2D semi-infinite domain in x and y directions.</figcaption></figure>

### A.4. 2D zero flux, finite domain in x, finite or semi-infinite in y

This subsection describes analytical solutions for problems with different configurations (see Figures 26, 27). The RWPT algorithm proposed further below must be able to handle particles displacements in such configurations.

The solution of the initial-value problem finite in x and semi-infinite in y (see Figure 26) is the solution of the 1D finite in x problem multiplied by the solution of the 1D semi-infinite in y problem.

$$ \forall t \ge 0; \forall x_{1:2} \le x \le x_{2:3}; $$

$$ \frac{\theta}{M_0} C(x, y, t) = [G(y_0, 2Dt, y) + G(2y_{BC} - y_0, 2Dt, y)] \frac{\theta^x}{M_0} C_{BC12}^x(x, t). \quad (44) $$

The solution of the initial-value problem finite in x and y (see Figure 27) is the solution of the 1D finite in x multiplied by the solution of the 1D finite in y problem.
<figure><img src="image_30.png" /><figcaption>Figure 26. 2D finite domain in x and semi-infinite in y.</figcaption></figure>

<figure><img src="image_31.png" /><figcaption>Figure 27. 2D finite domain in x and y.</figcaption></figure>

$$ \forall t \ge 0; \forall x_{1:2} \le x \le x_{2:3}; $$

$$ \frac{\theta}{M_0} C(x, y, t) = \frac{\theta^y}{M_0} C_{BC12}^y(y, t) \frac{\theta^x}{M_0} C_{BC12}^x(x, t). \quad (45) $$

References

[1] S. Whitaker, *The Method of Volume Averaging, Theory and Applications of Transport in Porous Media*, vol. 13, Springer, Dordrecht, 1999.

[2] M. Quintard, S. Whitaker, “One- and two-equation models for transient diffusion processes in two-phase systems”, in *Advances in Heat Transfer*, vol. 23, Academic Press, 1993, p. 369-464.

[3] E. Sanchez-Palencia, “Homogenization method for the study of composite media”, in *Asymptotic Analysis II*, Springer, Berlin, Heidelberg, 1983, p. 192-214.
[4] C. Moyne, S. Didierjean, H. Amaral Souto, O. da Silveira, "Thermal dispersion in porous media: one-equation model", *Int. J. Heat Mass Transfer* **43** (2000), p. 3853-3867.

[5] H. Brenner, "Dispersion resulting from flow through spatially periodic porous media", *Phil. Trans. R. Soc. Lond. A* **297** (1980), no. 1430, p. 81-133.

[6] J. L. Auriault, P. M. Adler, "Taylor dispersion in porous media: analysis by multiple scale expansions", *Adv. Water Resour.* **18** (1995), no. 4, p. 217-226.

[7] W. E. B. Engquist, X. Li, W. Ren, E. Vanden-Eijnden, "Heterogeneous multiscale methods: a review", *Commun. Comput. Phys.* **2** (2007), no. 3, p. 367-450.

[8] M. Vasilyeva, E. T. Chung, Y. Efendiev, A. Tyrylgin, "A three-level multi-continua upscaling method for flow problems in fractured porous media", *Commun. Comput. Phys.* **27** (2020), no. 2, p. 619-638.

[9] B. P. Muljadi, B. Bijeljic, M. J. Blunt, A. Colbourne, A. J. Sederman, M. D. Mantle, L. F. Gladden, "Modelling and upscaling of transport in carbonates during dissolution: validation and calibration with NMR experiments", *J. Contam. Hydrol.* **212** (2018), p. 85-95.

[10] A. J. Guttmann, T. Kennedy, "Self-avoiding walks in a rectangle", *J. Eng. Math.* **84** (2014), p. 201-208.

[11] H. Oukili, R. Ababou, G. Debenest, B. Noetinger, "Random Walks with negative particles for discontinuous diffusion and porosity", *J. Comput. Phys.* **396** (2019), p. 687-701.

[12] B. Noetinger, D. Roubinet, A. Russian, T. Le Borgne, F. Delay, M. Dentz et al., "Random walk methods for modeling hydrodynamic transport in porous and fractured media from pore to reservoir scale", *Transp. Porous Media* **115** (2016), no. 2, p. 345-385.

[13] B. Noetinger, T. Estebenet, "Up-scaling of double porosity fractured media using continuous-time random walks methods", *Transp. Porous Media* **39** (2000), no. 3, p. 315-337.

[14] B. Noetinger, T. Estebenet, M. Quintard, "Up-scaling flow in fractured media: equivalence between the large scale averaging theory and the continuous time random walk method", *Transp. Porous Media* **43** (2001), no. 3, p. 581-596.

[15] B. Noetinger, T. Estebenet, P. Landereau, "A direct determination of the transient exchange term of fractured media using a continuous time random walk method", *Transp. Porous Media* **44** (2001), no. 3, p. 539-557.

[16] R. Mauri, "Heat and mass transport in random velocity fields with application to dispersion in porous media", *J. Eng. Math.* **29** (1995), p. 77-89.

[17] L. C. Evans, *Partial Differential Equations*, 2nd ed., American Mathematical Society, Providence, RI, 2010.

[18] H. S. Carslaw, J. C. Jaeger, *Conduction of Heat in Solids*, Clarendon, Oxford, 1959.

[19] G. J. M. Uffink, "A random-walk method for the simulation of macrodispersion in a stratified aquifer", in *Relation of Groundwater Quality and Quantity*, IAHS Publ., vol. 146, Int. Assoc. of Hydrol. Sci., Gentbrugge, Belgium, 1985, p. 103-114.

[20] E. M. LaBolle, G. E. Fogg, A. F. B. Tompson, "Random-walk simulation of transport in heterogeneous porous media: local mass-conservation problem and implementation methods", *Water Resour. Res.* **32** (1996), no. 3, p. 583-593.

[21] P. Salamon, D. Fernandez-Garcia, J. J. Gomez-Hernandez, "A review and numerical assessment of the random walk particle tracking method", *J. Contam. Hydrol.* **87** (2006), no. 3-4, p. 277-305.

[22] M. Bechtold, J. Vanderborght, O. Ippisch, H. Vereecken, "Efficient random walk particle tracking algorithm for advective-dispersive transport in media with discontinuous dispersion coefficients and water contents", *Water Resour. Res.* **47** (2011), article no. W10526.

[23] Y. Yong, X. Lou, S. Li, C. Yang, X. Yin, "Direct simulation of the influence of the pore structure on the diffusion process in porous media", *Comput. Math. Appl.* **67** (2014), no. 2, p. 412-423.

[24] J. A. Currie, "Gaseous diffusion in porous media: II. Dry granular materials", *Br. J. Appl. Phys.* **11** (1960), p. 318-324.

[25] R. J. Millington, J. P. Quirk, "Permeability of porous solids", *Trans. Faraday Soc.* **57** (1961), p. 1200-1207.

[26] D. Ryan, R. G. Carbonell, S. Whitaker, "A theory of diffusion and reaction in porous media", in *AIChE Symposium Series, No. 202* (P. Sfroeve, W. J. Ward, eds.), vol. 77, 1981, p. 46-62.

[27] J. C. Maxwell, *Treatise on Electricity and Magnetism*, 2nd ed., vol. I, Clarendon Press, Oxford, 1881.

[28] H. L. Weissberg, J. C. Maxwell, *Treatise on Electricity and Magnetism*, 2nd ed., vol. I, Clarendon Press, Oxford, 1881, "Effective diffusion coefficients in porous media", *J. Appl. Phys.* **34** 1963, p. 2636-2639.

[29] N. Wakao, J. M. Smith, "Diffusion in catalyst pellets", *Chern. Eng. Sci.* **17** (1962), p. 825-834.

[30] J. H. Kim, J. A. Ochoa, S. Whitaker, "Diffusion in anisotropic porous media", *Transp. Porous Media* **2** (1987), no. 4, p. 327-356.

[31] A. G. Hunt, M. Sahimi, "Flow, transport, and reaction in porous media: percolation scaling, critical-path analysis, and effective medium approximation", *Rev. Geophys.* **55** (2017), p. 993-1078.