# Modeling diffusion in discontinuous media under generalized interface conditions: theory and algorithms

Elisa Baioni, Antoine Lejay, Géraldine Pichot, Giovanni Michele Porta

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

## To cite this version:

Elisa Baioni, Antoine Lejay, Géraldine Pichot, Giovanni Michele Porta. Modeling diffusion in discontinuous media under generalized interface conditions: theory and algorithms. 2023. hal-04166559v2

HAL Id: hal-04166559

https://inria.hal.science/hal-04166559v2

Preprint submitted on 5 Aug 2023

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
Modeling diffusion in discontinuous media under
generalized interface conditions: theory and
algorithms

Elisa Baioni¹, Antoine Lejay², Géraldine Pichot³, and Giovanni Michele
Porta¹

¹Politecnico di Milano, Department of Civil and Environmental Engineering, 20133 Milan (IT)

²Université de Lorraine, CNRS, IECL, Inria, F-5400 Nancy, France; ORCID: 0000-0003-0406-9550

³Inria - 2 rue Simone Iff, 75589 Paris, France and Université Paris-Est, CERMICS (ENPC), 6 et 8 av. Blaise Pascal, 77455 Marne-la-Vallée Cedex 2, France

Abstract

Diffusive transport in media with discontinuous properties is a challenging problem that arises in many applications. This paper presents a novel analytical expression from the method of images to simulate diffusive processes, such as mass or thermal transport, in discontinuous media with generalized boundary conditions at the discontinuity interface relying upon the random walk method. The analytical expressions are used to formulate a generalization of the existing Skew Brownian Motion, HYMLA and Uffink's method, here named as GSBM, GHYMLA and GUM respectively, to handle generic interface conditions. The algorithms are tested by simulating transport in a bimaterial medium with piece-wise constant properties. The results indicate that the GUM algorithm provides the best performance in terms of accuracy and computational cost. The methods proposed can be applied for simulation of a wide range of differential problems.

**Keywords–** Discontinuous media, diffusive transport, analytical solution, numerical simulation, random walk methods

# 1 Introduction

Modeling diffusion processes in discontinuous media is a challenging problem of both active
experimental [6] and theoretical interest [26, 1, 4, 14, 16, 12, 22]. Diffusive discontinuities
commonly occur at the variation of the physical properties of the medium (composite or
fractured media, geological repositories, layered porous media) or changes in phase, such
as the interface between surface water and groundwater. Processes involving pure diffusion
include heat diffusion in material (Fourier's law) and pressure diffusion in fluid-filled porous
media.

Random-walk particle-tracking (RWPT) methods are valid numerical approaches to deal
with diffusive discontinuity interfaces. The RWPT methods have been used to model the
transport of mass, heat and other physical quantities in solids, fluids and porous media using
particles as units of transported quantity (solute concentration, fluid pressure, heat, thermal
energy). This approach avoids some problems associated with the Eulerian methods relying
upon partial differential equations, such as mass balance error and numerical instability, and
can be computationally effective since each particle's behavior is independent of the others.
As demonstrated in [13], under diffusive discontinuities the standard random walk does not
guarantee mass balance and does not simulate the exact spatial distribution of particles. A
specific study for predicting the correct particle's behavior close to the discontinuity interface
is necessary. Two main conceptual schemes are used to predict the particle's displacement
in this context: (i) the interpolation method [13, 23, 11], (ii) the partial reflection [12, 26,
18, 1, 5, 19, 20, 3] approach, including a nonlinearly decomposed step [4]. The first class
smooths out the discontinuities by interpolating the discontinuous values on the numerical grid.
This method requires a high spatial resolution grid close to the discontinuity and a small time
step, which results in a significant increase in the computational cost. The partial reflection
approach was introduced by Uffink [26] and considers the discontinuity as a semi-reflective
barrier: the particles can cross the interface or be reflected based on the probability depending
on the diffusive properties of the discontinuous medium. This approach does not require grid
refinement and is less computationally expensive than the interpolation method. The partial
reflection method has been largely applied to simulate diffusive transport in heterogeneous
media [12, 27, 17, 24, 31]. In this regard, the main computational schemes are recalled in
the following. Uffink [26] suggested that in a bilayered medium, where a constant diffusion
coefficient characterizes each layer, a fraction of the particles moving from one layer to another
must be reflected. Cordes et al. [9] revisited Uffink's method by retaining the idea of the
semi-reflective barrier but modifying the fraction of reflected particles. Ackerer introduced
the “splitting up jump” approach in which the particle jump is divided into two components:
(1) the displacement occurring within one layer until the discontinuity interface, and (2) the
motion from discontinuity to the other layer. Hoteit, Mose, Younes, Lehmann, and Ackerer
[12] presented a new reflection method, called HMYLA, which combines the jump splitting
approach formulated by Ackerer with Semra's scheme [25] for the computation of the particle's
behavior at the discontinuity. Lastly, the Skew Brownian motion (SBM) model implemented
by Lejay and Pichot [16, 17] describes the particle motion around the discontinuity interface
by a stochastic process grounded on the Skew Brownian motion. According to [17], the SBM
is the algorithm providing higher accuracy in simulating diffusive transport in a discontinuous
system.

Mass-transfer particle tracking (MTPT) methods have recently been combined with partial
reflection methods to solve diffusive transport in discontinuous diffusion media. These
Lagrangian approaches simulate diffusion by both random walk and diffusive mass transfer
and differ from conventional RWPT methods since particle masses can be transferred among
particles. Schmidt et al. [24] implemented an MT algorithm that incorporates a semi-analytical
solution to deal with diffusive discontinuities and that can be straightforwardly extended to
higher dimensions and complicated subdomain interfaces. Lastly, the restriction on a small
time step in the reflection method, needed to ensure correct particle behavior close to the
discontinuity, is overcome in [21] through the implementation of a new “Stop&Go” (i.e.,
two steps) RWPT method combining partial reflection schemes, as in [16, 12], with adaptive
negative mass particles.

The approaches mentioned above solve diffusive problems in media characterized by a dis-
continuous diffusivity with a continuous solution imposed at the discontinuity interface, a
condition required to model mass transport problems in heterogeneous media (e.g., porous
media). Other physical problems, such as thermal diffusion, require different interface condi-
tions. Here, we propose a novel analytical expression from the method of the images [26] for
simulating transport processes, such as mass or thermal diffusion, in discontinuous media with
generalized boundary conditions at the discontinuity relying upon the random walk method.
Different boundary conditions at the discontinuity are obtained by computing four constants
depending on the physical properties of the medium, the latter being assumed piece-wise
constant. Based on these analytical results we formulate a generalized version of three selected
algorithms, i.e., the Uffink's, SBM and the HYMLA algorithm, thus deriving generalized
Uffink method (GUM), SBM (GSBM) and HYMLA (GHYMLA) approaches.

The paper is structured as: Section 2 introduces an analytical formulation from the method of the images to simulate a one-dimensional diffusion process in a bimaterial medium with piece-wise constant diffusivity and a discontinuity at 0. The approach is presented for the point source located in both the negative (Section 2.1) and positive (Section 2.2) region of the domain. The GUM method is illustrated in Section 3 while the GSBM and GHYMLA are described in Section 3.2 and in Appendix C, along with their corresponding computational algorithms. These methods solve diffusive problems in discontinuous media with generalized interface conditions. The benchmark test used to evaluate the algorithms' performance is described in Section 4, as well as a discussion of the numerical results. The conclusions drawn from our study are reported in Section 5.
## 2 Analytic expression from the method of images

We present our methodology considering the simulation through a random walk of a one-dimensional diffusion process in a bimaterial medium with a piece-wise constant diffusivity $D$ discontinuous at $X_I = 0$. The discontinuity interface is perceived as a semi-reflective barrier and divides the domain in $R_1$ for $x < 0$ and $R_2$ for $x > 0$ with diffusivity equal to $D_1$ and $D_2$, respectively.

More precisely, we consider solving the family of problems

$$ \partial_t q(t, x, y) = \nabla_y(D(y)\nabla_y q(t, x, y)), \quad (1) $$

$$ \nu_1 q(t, x, 0^{-}) = \nu_2 q(t, x, 0^{+}), \quad (2) $$

$$ \kappa_1 \nabla_y q(t, x, 0^{-}) = \kappa_2 \nabla_y q(t, x, 0^{+}), \quad (3) $$

$$ q(t, x, y) \underset{t \to 0}{\longrightarrow} \delta_y(x), \quad (4) $$

with

$$ D(y) := \begin{cases} D_1 & y < 0, \\ D_2 & y > 0. \end{cases} $$

The interface parameters $\nu_i, \kappa_i$, with $i = 1, 2$, are constants characterizing each layer, such as physical properties of the medium. The function $(t, x, y) \mapsto q(t, x, y)$ of solutions to (1)-(4) is called a fundamental solution.

The fundamental solution $q$ is derived from the method of images [10] for the medium defined here. Available algorithms such as the Uffink method (UM) [27], the SBM [16, 17] and the HYMLA [12] can be derived from the method of images under the condition where $\nu_1 = \nu_2$ and $\kappa_i = D_i$. We extend and generalize these methods by relaxing these assumptions and thus considering (2)-(3).

The method of images consists in constructing the fundamental solution as a superposition of simpler fundamental solutions with different source points. For this, we introduce the Gaussian densities

$$ p_i(t, x, y) := \frac{1}{\sqrt{4\pi D_i t}} \exp\left(\frac{-|y - x|^2}{4D_i t}\right) \quad (5) $$

with $i = 1, 2$. For each $x, (t, y) \mapsto p_i(t, x, y)$ solves the heat equation

$$ \partial_t p(t, x, y) = D_i \Delta_y p_i(t, x, y) \quad \text{with} \quad p_i(t, x, y) \underset{t \to 0}{\longrightarrow} \delta_x(y). \quad (6) $$
## 2.1 Method of images for a negative source

We now consider a point source located in region $R_1$, i.e., at $x < 0$, here termed as *negative source*. Following [16], our guess for the fundamental solution $q$ is

$$ q(t, x, y) := \begin{cases} p_1(t, x, y) + Ap_1(t, -x, y) & \text{if } y < 0, \\ Bp_2(t, \beta x, y) & \text{if } y \ge 0 \end{cases} \qquad (7) $$

with $\beta := \frac{\sqrt{D_2}}{\sqrt{D_1}}$ and $A$ and $B$ fixed and corresponding to the fraction of particles that are reflected and cross the discontinuity interface, respectively. Since $x < 0$, we obtain that $q$ satisfies (1)-(4), where $\delta_y(x)$ is the delta function of the particles concentration in $x$.

The parabolic problem (1)-(4) is not properly defined unless one specifies the constants $A$ and $B$, which depend on the behavior of the solution at the interface. At the discontinuity interface $X_I = 0$ we find

$$ \left\{ \begin{aligned} q(t, x, 0^{-}) &= \frac{1+A}{\sqrt{4\pi D_1 t}} \exp\left(\frac{-x^2}{4D_1 t}\right), \\ q(t, x, 0^{+}) &= \frac{B}{\sqrt{4\pi D_2 t}} \exp\left(\frac{-\beta^2 x^2}{4D_2 t}\right) = \frac{B}{\beta\sqrt{4\pi D_1 t}} \exp\left(\frac{-x^2}{4D_1 t}\right), \end{aligned} \right. \qquad (8) $$

and

$$ \left\{ \begin{aligned} \nabla_y q(t, x, 0^{-}) &= \frac{(1-A)x}{D_1 t \sqrt{4\pi D_1 t}} \exp\left(\frac{-x^2}{4D_1 t}\right), \\ \nabla_y q(t, x, 0^{+}) &= \frac{B\beta x}{D_2 t \sqrt{4\pi D_2 t}} \exp\left(\frac{-\beta^2 x^2}{4D_2 t}\right) = \frac{Bx}{\beta^2 D_1 t \sqrt{4\pi D_1 t}} \exp\left(\frac{-x^2}{4D_1 t}\right). \end{aligned} \right. \qquad (9) $$

We now compute $A$ and $B$ to consider the interface conditions in Equations (2)-(3). According to Appendix A, with

$$ \mu := \frac{\nu_2 - \nu_1}{\nu_2 + \nu_1} \quad \text{and} \quad \gamma := \frac{\kappa_2 - \kappa_1}{\kappa_2 + \kappa_1}, \qquad (10) $$

we may rewrite these conditions as

$$ (1 - \mu)q(t, x, 0^{-}) = (1 + \mu)q(t, x, 0^{+}), \qquad (11) $$

$$ (1 - \gamma)\nabla_y q(t, x, 0^{-}) = (1 + \gamma)\nabla_y q(t, x, 0^{+}) \qquad (12) $$

so that only two parameters $\mu$ and $\gamma$ are involved.
After some eliminations in Equations (8)-(9), Equations (11)-(12) reduce to

$$ (1 - \gamma)\beta^2(1 - A) = (1 + \gamma)B, \qquad (13) $$

$$ (1 - \mu)\beta(1 + A) = (1 + \mu)B. \qquad (14) $$

Therefore, the constants A and B, with $A \in (-1, 1)$ and $B > 0$, write as

$$ A = \frac{\Gamma - 1}{\Gamma + 1}, \qquad (15) $$

$$ B = \frac{1 - \gamma}{1 + \gamma} \frac{2\beta^2}{1 + \Gamma}, \qquad (16) $$

where

$$ \Gamma := \beta \frac{1 - \gamma}{1 + \gamma} \frac{1 + \mu}{1 - \mu} = \frac{\sqrt{D_2} \nu_2 \kappa_1}{\sqrt{D_1} \nu_1 \kappa_2}. $$

## 2.2 Method of images for a positive source

We now consider a positive source $x > 0$ and follow a similar procedure as the one described in Section 2.1. The fundamental solution $q$ is defined as

$$ q(t, x, y) := \begin{cases} \tilde{B}p_1(t, \tilde{\beta}x, y) & \text{if } y < 0, \\ p_2(t, x, y) + \tilde{A}p_2(t, -x, y) & \text{if } y \ge 0 \end{cases} \qquad (17) $$

with $\tilde{\beta} := \frac{\sqrt{D_1}}{\sqrt{D_2}} = \frac{1}{\beta}$ and $\tilde{A}, \tilde{B}$ fixed and corresponding to the fraction of particles that are reflected and cross the discontinuity interface, respectively. At the discontinuity interface we find

$$ \left\{ \begin{aligned} q(t, x, 0^{-}) &= \frac{\tilde{B}}{\sqrt{4\pi D_1 t}} \exp\left(\frac{-\tilde{\beta}^2 x^2}{4D_1 t}\right) = \frac{\tilde{B}}{\tilde{\beta}\sqrt{4\pi D_2 t}} \exp\left(\frac{-x^2}{4D_2 t}\right), \\ q(t, x, 0^{+}) &= \frac{1 + \tilde{A}}{\sqrt{4\pi D_2 t}} \exp\left(\frac{-x^2}{4D_2 t}\right) \end{aligned} \right. \qquad (18) $$

and

$$ \left\{ \begin{aligned} \nabla_y q(t, x, 0^{-}) &= \frac{\tilde{B}x}{D_1 t \sqrt{4\pi D_1 t}} \exp\left(\frac{-\tilde{\beta}^2 x^2}{4D_1 t}\right) = \frac{\tilde{B}x}{\tilde{\beta}^2 D_2 t \sqrt{4\pi D_2 t}} \exp\left(\frac{-x^2}{4D_2 t}\right), \\ \nabla_y q(t, x, 0^{+}) &= \frac{(1 - \tilde{A})x}{D_2 t \sqrt{4\pi D_2 t}} \exp\left(\frac{-x^2}{4D_2 t}\right). \end{aligned} \right. \qquad (19) $$
Using $\mu$ and $\gamma$ defined in (10) and proceeding as for negative sources,

$$ (1 - \mu)\tilde{B} = (1 + \mu)(1 + \tilde{A})\tilde{\beta}, \qquad (20) $$

$$ (1 - \gamma)\tilde{B} = (1 + \gamma)(1 - \tilde{A})\tilde{\beta}^2. \qquad (21) $$

In summary, for positive source in $x > 0$ the constants $\tilde{A}$ and $\tilde{B}$, with $\tilde{A} \in (-1, 1)$ and $B > 0$, are formulated as

$$ \tilde{A} = \frac{\tilde{\Gamma} - 1}{\tilde{\Gamma} + 1} = \frac{1 - \Gamma}{1 + \Gamma} = -A, \qquad (22) $$

$$ \tilde{B} = \frac{1 + \gamma}{1 - \gamma} \frac{2\tilde{\beta}^2}{1 + \tilde{\Gamma}} = \frac{\sqrt{D_1} D_1 \nu_2 \kappa_2}{\sqrt{D_2} D_2 \nu_1 \kappa_1} B = \frac{1}{B} \frac{4\Gamma}{(1 + \Gamma)^2} \qquad (23) $$

with

$$ \tilde{\Gamma} := \tilde{\beta} \frac{1 + \gamma}{1 - \gamma} \frac{1 - \mu}{1 + \mu} = \frac{1}{\Gamma}. $$

## 2.3 Mass conservation

We derived in (7) and (17) an analytic expression for the solution $q$ to (1)-(4) with the interface conditions (2)-(3). This problem involves the 6-uple ($D_1, D_2, \nu_1, \nu_2, \kappa_1, \kappa_2$) of parameters, but depends actually only on the 4-uple ($D_1, D_2, \nu_1/\nu_2, \kappa_1/\kappa_2$), or, using (10), ($D_1, D_2, \mu, \gamma$). The mass conservation means that

$$ \int_{-\infty}^{+\infty} q(t, x, y) \, dy = 1 \qquad (24) $$

for any source point $x$ and any $t \ge 0$. Such a property is essential to interpret $q(t, x, \cdot)$ as the density of a random variable, and then use this to develop random walk schemes.

As proven in Appendix B, for $x < 0$, the condition $A + B = 1$ is equivalent to the mass conservation. This condition is met if and only if

$$ \beta^2 = \frac{D_2}{D_1} = \frac{1 + \gamma}{1 - \gamma} = \frac{\kappa_2}{\kappa_1}. \qquad (25) $$

The same result holds for both positive and negative source solutions. Note that mass conservation depends only on conditions on $\kappa_i$ and $D_i$, $i = 1, 2$, but is independent of $\nu_i$.
## 2.4 The Kolmogorov backward equation

The fundamental solution $q$ was constructed as solving a family of partial differential equations with respect to the variable $y$. Eq. (1)-(4) is referred as *Kolmogorov forward equation*, the *Master equation* or the *Fokker-Planck equation*.

The same problem can be reformulated writing $q$ as a function of $x$, deriving the so-called *Kolmogorov backward equation*, where the terms *forward* and *backward* refers to the variables $y$ and $x$ on which the operators act. This derivation yields:

$$ \partial_t q(t, x, y) = \nabla_x (D(x) \nabla_x q(t, x, y)), \quad (26) $$

$$ \frac{D_1}{\kappa_1} q(t, 0^-, y) = \frac{D_2}{\kappa_2} q(t, 0^+, y), \quad (27) $$

$$ \frac{D_1}{\nu_1} \nabla_y q(t, 0^-, y) = \frac{D_2}{\nu_2} \nabla_y q(t, 0^+, y), \quad (28) $$

$$ q(t, x, y) \xrightarrow{t \to 0} \delta_y(x). \quad (29) $$

In this formulation the interface conditions are different with respect to their forward formulation (2)-(3). Although establishing (26)-(29) can be performed from the expressions (7) and (17), this obeys to a general principle [30] by identifying the dual operator of $\nabla_y(D(y)\nabla_y\cdot)$ with the interface conditions (2)-(3). We refer to [7] for the derivation of our interface conditions.

We remark that the Condition (25) which ensures the mass conservation is equivalent to the continuity of $x \mapsto q(t, x, y)$ at 0 in (27).

Under Condition (25), when interpreting $q(t, x, \cdot)$ as a probability density of a particle position initially at $x$ at time $t$, then we may set $x = 0$. The quantities

$$ \mathbb{P}_+ := \int_0^L q(t, 0, y) \, dy = \frac{1 + \tilde{A}}{2} = \frac{1 - A}{2} = \frac{B}{2} \quad \text{and} \quad \mathbb{P}_- := \int_{-L}^0 q(t, 0, y) \, dy = 1 - \frac{B}{2}. $$

represent the probabilities for the particle starting from 0 to displace on the right and left side of the interface, respectively.

The mass conservation condition is $A + B = 1$. Since $A \in (-1, 1)$, $B \in (0, 2)$. We then define

$$ \theta_g := B - 1 = -A \in (-1, 1) \text{ so that } \mathbb{P}_+ = \frac{1 + \theta_g}{2} \text{ and } \mathbb{P}_- = \frac{1 - \theta_g}{2}. \quad (30) $$

In this case, when the mass is conserved,

$$ \mathbb{P}_+ = \frac{B}{2} = \frac{1}{1 + \Gamma}, \quad \mathbb{P}_- = \frac{\Gamma}{1 + \Gamma}, \quad B = 1 + \theta_g \text{ and } A = -\theta_g. $$

This result implies that $\theta_g$ provides a convenient reparameterization of the pair $(A, B)$ with
the constraint $A + B = 1$.

# 3 Approximation algorithms for generalized interface conditions

Since $q(t, x, y) > 0$ and considering mass conservation, $q(t, x, \cdot)$ can be interpreted as a probability. Adopting a Lagrangian perspective, we then sample particles' positions according to this. In the following sections we provide three algorithms to deal with particle displacements under the considered interface conditions. Note that in Algorithms 1, 2, 3 below, we denote by $N(0, 1)$ for the Gaussian distribution, $U(0, 1)$ for the Uniform distribution in $(0, 1)$ and $IG(\alpha, \beta)$ for the Inverse Gaussian distribution of parameters $(\alpha, \beta)$.

## 3.1 Generalized Uffink Method

An algorithm for the generalized Uffink method (GUM) can be developed from the method of images described in Sections 2.1-2.2. Here, we demonstrate how to approximate the probability transition density for both a negative and positive source.

The algorithm is then given in Algorithm 1.

### 3.1.1 Case of negative source

Let us consider a negative source, $x < 0$. The density $q$ writes as Equation (7) with $A, B$ defined according to Equations (15)-(16), respectively. To be consistent with the original Uffink method [27], the Gaussian density $p_1(\Delta t, x, y)$ is replaced by a uniform density $r_1(\Delta t, x, y)$ where $r_1(\Delta t, x, \cdot)$ is piecewise constant. Thus we obtain

$$r_1(\Delta t, x, y) = \frac{1}{2H_1} \mathbb{1}_{[x-H_1, x+H_1]}(y) \quad (31)$$

with $H_1 = \sqrt{6D_1\Delta t}$.

There are three possible cases to consider:

1. For $x + H_1 < 0$ the particle is far enough to not be affected by the diffusive discontinuity. Therefore, a uniformly distributed step is employed;

2. For $x, y < 0$ only the particles on the region $R_1(x < 0)$ are considered. The density $p_1(\Delta t, -x, y) = p_1(\Delta t, x, -y)$ is replaced by $r_1(\Delta t, x, -y) \mathbb{1}_{y \le 0}$;
3. For $x < 0$, $y > 0$ the particle crosses the interface during the time step. The density $p_2(\Delta t, \beta x, y)$ is replaced by $\frac{1}{\beta}r_1(\Delta t, x, y/\beta)1_{y \ge 0}$ since

$$p_2(\Delta t, \beta x, y) = \frac{1}{\sqrt{4\pi D_2 t}} \exp\left(\frac{-(y/\beta - x)^2 \beta^2}{4D_2 t}\right) = \frac{1}{\beta} p_1\left(\Delta t, x, \frac{y}{\beta}\right). \quad (32)$$

In summary, for a negative source the density $q(\Delta t, x, y)$ can be approximated by $r(\Delta t, x, y)$ defined as

$$r(\Delta t, x, y) := \begin{cases} r_1(\Delta t, x, y) + A r_1(\Delta t, x, -y) & \text{if } y < 0, \\ \frac{B}{\beta} r_1(\Delta t, x, y/\beta) & \text{if } y \ge 0. \end{cases} \quad (33)$$

Since $p_1, p_2$ are replaced by $r_1, r_2$ and the mass is conserved, $y \mapsto r(\Delta t, x, y)$ is a density of a random variable. A random variable from the density $r(\Delta t, x, \cdot)$ can be simulated by considering the splitting into three regions $J_1 = [x - H_1, -(x + H_1)]$, $J_2 = [-(x + H_1), 0]$, $J_3 = [0, (x + H_1)\beta]$, and drawing a uniform variable with step-wise constant probability on each of these zones as displayed in Figure 1. If $V$ has a density $r(\Delta t, x, \cdot)$, then

$$\mathbb{P}_1 = \mathbb{P}[V \in J_1] = \frac{-2x}{2H_1}, \quad (34)$$

$$\mathbb{P}_2 = \mathbb{P}[V \in J_2] = (x + H_1) \frac{1 + A}{2H_1}, \quad (35)$$

$$\mathbb{P}_3 = \mathbb{P}[V \in J_3] = (x + H_1) \frac{B}{2H_1}. \quad (36)$$

<figure><img src="image_3.png" /><figcaption>Figure 1: Probability density of $r(\Delta t, x, \cdot)$ for a negative source.</figcaption></figure>

The distribution of $V$ given $V \in (J_1|J_2|J_3)$ is uniform over $J_1, J_2, J_3$ and can be drawn according to $r(\Delta t, x, \cdot)$ from the knowledge of $A, B$ (see Equations (15)-(16)), the particle's position $x$ at the time step $\Delta t$ and $\beta$, the latter being independent of $x$.

We see from Figure 1 that if $(1 + A)\beta = B$, then we could merge $J_2$ and $J_3$ in a single step.
This corresponds to the original Uffink algorithm [27].

### 3.1.2 Case of positive source

We now demonstrate how to approximate the density $q(\Delta t, x, \cdot)$ (see Equation (17)) for source in $x > 0$. Similarly to Section 3.1.1, we may replace the Gaussian density $p_2(\Delta t, x, y)$ with the uniform density $\tilde{r}_2(\Delta t, x, y)$ as in Equation (37), where $\tilde{r}_2(\Delta t, x, \cdot)$ is constant.

$$ \tilde{r}_2(\Delta t, x, y) = \frac{1}{2H_2} \mathbb{1}_{[x-H_2, x+H_2]}(y) \quad (37) $$

with $H_2 = \sqrt{6D_2\Delta t}$. Let us consider three scenarios:

1. If $x - H_2 > 0$, we keep only the uniform step;

2. If $x, y > 0$, the density $p_2(\Delta t, -x, y) = p_2(\Delta t, x, -y)$ is replaced by $\tilde{r}_2(\Delta t, x, -y) \mathbb{1}_{y \ge 0}$ since only the particles in $R_2$ are taken into account;

3. If $x > 0, y < 0$, the density $p_1(\Delta t, \tilde{\beta}x, y)$ is substituted by $\frac{1}{\tilde{\beta}}\tilde{r}_2(\Delta t, x, y/\tilde{\beta}) \mathbb{1}_{y < 0}$ since

$$ p_1(\Delta t, \tilde{\beta}x, y) = \frac{1}{\sqrt{4\pi D_1 t}} \exp \left( -\frac{(y/\tilde{\beta} - x)^2 \tilde{\beta}^2}{4D_1 t} \right) = \frac{1}{\tilde{\beta}} p_2 \left( \Delta t, x, \frac{y}{\tilde{\beta}} \right). \quad (38) $$

Thus, the function $\tilde{r}(\Delta t, x, y)$ is defined as

$$ \tilde{r}(\Delta t, x, y) := \begin{cases} \frac{\tilde{B}}{\tilde{\beta}} \tilde{r}_2 \left( \Delta t, x, \frac{y}{\tilde{\beta}} \right) & \text{if } y < 0, \\ \tilde{r}_2(\Delta t, x, y) + \tilde{A} \tilde{r}_2(\Delta t, x, -y) & \text{if } y \ge 0. \end{cases} \quad (39) $$

As indicated in Figure 2 we can simulate a random variate from the density $\tilde{r}(\Delta t, x, \cdot)$ by considering three zones $\tilde{J}_1 = [(x - H_2)\tilde{\beta}, 0]$, $\tilde{J}_2 = [0, -(x - H_2)]$, $\tilde{J}_3 = [-(x - H_2), (x + H_2)]$ and by drawing a uniform variable with the correct probability for each region. Let $\tilde{V}$ have the density $\tilde{r}(\Delta t, x, \cdot)$, the probability of each zone is expressed as

$$ \mathbb{P}_1 = \mathbb{P}[\tilde{V} \in \tilde{J}_1] = -(x - H_2) \frac{\tilde{B}}{2H_2}, \quad (40) $$

$$ \mathbb{P}_2 = \mathbb{P}[\tilde{V} \in \tilde{J}_2] = -(x - H_2) \frac{1 + \tilde{A}}{2H_2}, \quad (41) $$

$$ \mathbb{P}_3 = \mathbb{P}[\tilde{V} \in \tilde{J}_3] = \frac{2x}{2H_2}. \quad (42) $$
<figure><img src="image_4.png" /><figcaption>Figure 2: Probability density of $\tilde{r}(\Delta t, x, \cdot)$ for a positive source.</figcaption></figure>

In conclusion, $\tilde{V}$ is drawn according to $\tilde{r}(\Delta t, x, \cdot)$ by knowing $\tilde{A}, \tilde{B}, H_2, \tilde{\beta}$, the particle's position $x$ and the time step $\Delta t$.

## 3.2 Generalized SBM (GSBM)

We call GSBM the SBM (Skew Brownian motion) method proposed in [16, 17] with a generalized parameter $\theta_g$ to deal with arbitrary interface conditions. Similarly, we propose a modification of the HYMLA algorithm introduced by Hoteit, Younes, Mose, Lehmann and Ackerer [12]. The generalized interface conditions (2)-(3) imposed at the discontinuity interface require a new formulation of the parameter $\theta$ employed in both the original SBM and HYMLA algorithms. The new parameter $\theta$, here referred as $\theta_g$, for the system described in Section 2 is demonstrated in the following [27]. Details on the GSBM and HYMLA algorithms are provided in Appendix C.

### 3.2.1 Density of the Generalized SBM

For some parameter $\theta_g \in (-1, 1)$, we introduce

$$p^{\dagger}(t, x, y) := \frac{\exp\left(-\frac{(x-y)^2}{2t}\right)}{\sqrt{2\pi t}} + \theta_g \operatorname{sgn}(y) \frac{\exp\left(-\frac{(|x|+|y|)^2}{2t}\right)}{\sqrt{2\pi t}}. \quad (43)$$

This density satisfies (1)-(4) with

$$D_1 = D_2 = \frac{1}{2}, \quad \frac{\nu_1}{\nu_2} = \frac{1 - \theta_g}{1 + \theta_g} \quad \text{and} \quad \kappa_1 = \kappa_2$$

or, using the parametrization with $\mu$ and $\gamma$ of (10),

$$D_1 = D_2 = \frac{1}{2}, \quad \mu = \theta_g \quad \text{and} \quad \gamma = 0.$$
<caption>**Algorithm 1:** GUM algorithm</caption>

**Input:** An initial position $X_t = x$ and a time step $\Delta t > 0$ with an interface at $X_I = 0$.

**Output:** A position $X_{t+\Delta t}$ according to the GUM algorithm.

Compute $H_1 = \sqrt{6D_1\Delta t}$ and $H_2 = \sqrt{6D_2\Delta t}$;

if $x + H_1 < 0$ then

/* The interface is not crossed: uniform step */

Generate a random variate $V \sim U(0, 1)$;

return $x + (2V - 1)H_1$

else

if $x < 0$ then

/* Negative source and the interface is crossed: compute $P_1$ (34), $P_2$ (35) and $P_3$ (36) */

Compute $x_L = x - H_1$, $x_M = -x - H_1$ and $x_R = (x + H_1) \beta$;

$P_1 = \frac{1}{2H_1}(x_M - x_L)$, $P_2 = \frac{1+A}{2H_1}(-x_M)$ and $P_3 = \frac{B}{2\beta H_1}x_R$;

Generate random variates $U \sim U(0, 1)$ and $V \sim U(0, 1)$;

if $U \le P_1$ then

return $x_L + (x_M - x_L)V$;

else if $P_1 < U \le P_1 + P_2$ then

return $x_M V$;

else

return $x_R V$;

else

if $x - H_2 > 0$ then

/* The interface is not crossed: uniform step */

Generate a random variate $V \sim U(0, 1)$;

return $x + (2V - 1)H_2$

else

/* Positive source and the interface is crossed: compute $P_1$ (40), $P_2$ (41) and $P_3$ (42) */

Compute $x_L = (x - H_2) \tilde{\beta}$, $x_M = -(x - H_2)$ and $x_R = x + H_2$;

$P_1 = \frac{\tilde{\beta}}{2\beta H_2}(-x_L)$, $P_2 = \frac{1+A}{2H_2}(x_M)$, $P_3 = \frac{1}{2H_2}(x_R - x_M)$;

Generate random variates $U \sim U(0, 1)$ and $V \sim U(0, 1)$;

if $U \le P_1$ then

return $x_L V$;

else if $P_1 < U \le P_1 + P_2$ then

return $x_M V$;

else

return $x_M + (x_R - x_M)V$;
For general values of $D_1$ and $D_2$, we define the GSBM probability transition density

$$q^{\dagger}(t, x, y) := \frac{1}{\sqrt{2D(y)}} p^{\dagger} \left( t, \frac{x}{\sqrt{2D(x)}}, \frac{y}{\sqrt{2D(y)}} \right) \quad (44)$$

with $D(x) := D_1$ if $x < 0$ and $D(x) := D_2$ if $x > 0$. Such a $q$ satisfies (1)-(4) as well with

$$\frac{\nu_1}{\nu_2} = \frac{1 - \theta_g}{1 + \theta_g} \quad \text{and} \quad \frac{\kappa_1}{\kappa_2} = \frac{D_1}{D_2}.$$

With such conditions, (25) is satisfied so that the mass is conserved.

Comparing (44) with (7) and (17), obtained from the method of images, $q(t, x, y) = q^{\dagger}(t, x, y)$ for any $t > 0$, $x, y \in \mathbb{R}$ if and only if

$$\theta_g = -A = \tilde{A}. \quad (45)$$

Besides the definition of this new parameter $\theta_g$, the algorithms of the GSBM and GHYMLA practically correspond to their original versions, as $p^{\dagger}(t, x, y)$ is the density of a Skew Brownian motion of parameter $\theta_g$ [15, 28].

# 4 Benchmark test

We test here the performance of the GUM, GSBM and GHYMLA method in replicating the analytical solution obtained in Section 2. We first introduce the simulation set up and then illustrate our numerical results.

## 4.1 Simulation set up

We design a Lagrangian method based on $N$ particles on the finite domain $[-L, L]$ with reflection boundary conditions (RBC) at $x = \pm L$ (i.e., no-flux boundary conditions). We employ the GUM Algorithm 1, GSBM Algorithm 2 and GHYMLA Algorithm 3 to simulate diffusive transfer in a one-dimensional bimaterial medium as in Figure 3. The domain consists of two regions divided by an interface $X_I = 0$, i.e., $R_1$ for $x < 0$ and $R_2$ for $x > 0$. The physical properties of the domain, which are embedded in the model parameters $\nu_1, \nu_2, \kappa_1, \kappa_2$ in (2)-(3), are constant within $R_1$ and $R_2$. The particles are initially located at $X_0 \in [-L, L]$. The aim of the test is to check if the particles are distributed according to the analytical $q$.

When we are restricted to a finite domain, the closed-form expression $q$ has to be replaced by a fundamental solution $q_L$ that still satisfies the Kolmogorov forward equation and the
<figure><img src="image_5.png" /><figcaption>Figure 3: Test case setup: the blue and red-colored regions indicate regions $R_1, R_2$, the orange point represents the point source location, RBC stands for reflective boundary conditions.</figcaption></figure>

Kolmogorov backward equation (26)-(29) with the additional conditions

$$ \nabla_x q_L(t, \pm L, y) = \nabla_y q_L(t, x, \pm L) = 0. $$

Moreover, in presence of boundaries at $\pm L$, $q(t, x, \cdot)$ gives an approximation of $q_L$ when $|x|^2/t \ll 1$ (for an interface at 0) and $t \ll 1$. The methods GUM, GSBM and GHYMLA allows one to sample particles' position according to the density. The corresponding constants are given in Table (1).

<table><caption>Table 1: Constants in the algorithms for generalized interface conditions</caption><thead><tr><th>Algorithm</th><th colspan="2">Parameters</th></tr></thead><tbody><tr><td rowspan="2">GUM</td><td>A = <sup>√κ<sub>1</sub>/ν<sub>1</sub> - √κ<sub>2</sub>/ν<sub>2</sub></sup>&frasl;<sub>√κ<sub>1</sub>/ν<sub>1</sub> + √κ<sub>2</sub>/ν<sub>2</sub></sub></td><td>B = <sup>2√κ<sub>2</sub>/ν<sub>2</sub></sup>&frasl;<sub>√κ<sub>1</sub>/ν<sub>1</sub> + √κ<sub>2</sub>/ν<sub>2</sub></sub></td></tr><tr><td>~A = <sup>√κ<sub>2</sub>/ν<sub>2</sub> - √κ<sub>1</sub>/ν<sub>1</sub></sup>&frasl;<sub>√κ<sub>1</sub>/ν<sub>1</sub> + √κ<sub>2</sub>/ν<sub>2</sub></sub></td><td>~B = <sup>2√κ<sub>1</sub>/ν<sub>1</sub></sup>&frasl;<sub>√κ<sub>1</sub>/ν<sub>1</sub> + √κ<sub>2</sub>/ν<sub>2</sub></sub></td></tr><tr><td>GSBM, GHYMLA</td><td colspan="2">θ<sub>g</sub> = <sup>√κ<sub>2</sub>/ν<sub>2</sub> - √κ<sub>1</sub>/ν<sub>1</sub></sup>&frasl;<sub>√κ<sub>1</sub>/ν<sub>1</sub> + √κ<sub>2</sub>/ν<sub>2</sub></sub></td></tr></tbody></table>

The position $X(0)$ of a particle is drawn with a given initial distribution $\nu$ over $[-L, L]$. Hence, position $X(t)$ is determined according to the transition probability density $q_L(t, X(0), \cdot)$. For this, thanks to the Chapman-Kolmogorov equation, we fix a time step $\Delta t = t/n$ and we draw iteratively $X((k + 1)\Delta t)$ from any suitable approximation of $q_L(t, X(k\Delta t), \cdot)$. Thus, $\{X(k\Delta t)\}_{k=0,...,n}$ performs a random walk, and the distribution of $X(t)$ is close to $q_L(t, X(0), \cdot)$. The computational domain comprises two main regions:

1. *interface layer*: it is an interval $I$, which we may assume as $I = [-\sqrt{6D_1\Delta t}, \sqrt{6D_2\Delta t}]$, containing the interface point $X_I = 0$ such that for each $x \notin I$, the probability that the particle crosses $X_I$ is below a given threshold. Here, the next position is sampled by using $q(\Delta t, x, \cdot)$, or any of its approximations.
2. boundary layer: is an interval $B_z$, typically chosen as $B_L = [L - \sqrt{6D_2\Delta t}, L]$, containing either $z = L$ or $z = -L$ such that for any $x \notin B$, the probability that the particle crosses $X_I$ is below a given threshold. In the boundary layer $B_L$, we sample the next position by performing a Gaussian step (when using GSBM) or a uniform step (when using GUM or GHYMLA) and applying a reflection around $L$ if the next position is outside the domain.

The size of the interfaces and boundary layers are ruled by $\Delta t$. We select $\Delta t$ so that $I \cap B_L = I \cap B_{-L} = \emptyset$. Outside these intervals, we perform an usual step using a Gaussian (when using GSBM) and or a uniform (when using GUM or GHYMLA) step.

For counting the particles around a given point, we decompose the domain $[-L, L]$ into $N_a$ non-overlapping cells $A_i := [x_{i-1}, x_i]$ for $i = 1, \dots, N_a$ where the interface $0 = x_i$ for some index $i$, $x_0 = -L$, $x_{N_a} = L$. We write $a_i := (x_{i-1} + x_i)/2$, the midpoint of the cell $A_i$, and $|A_i| := x_i - x_{i-1}$, the length of the cell $A_i$.

For each simulation, the average computational time $t_{cm}$ [s] is computed as the time to complete the entire simulation divided by the number of steps. Each simulation is performed in Matlab with an Intel(R) Xeon(R) CPU E5-2640 v3 @ 2.60GHz processor. The computational parameters and the size of the domain are reported in Table 2 and Table 3. All the quantities are assumed non-dimensional.

<table><caption>Table 2: Size of the domain and model parameters used in the tests.</caption><thead><tr><th></th><th>R<sub>1</sub></th><th>R<sub>2</sub></th><th>parameter</th></tr></thead><tbody><tr><td>ν</td><td>2.39 × 10<sup>-7</sup></td><td>8.30 × 10<sup>-7</sup></td><td>model parameter</td></tr><tr><td>κ</td><td>1.44 × 10<sup>-7</sup></td><td>2.78 × 10<sup>-7</sup></td><td>model parameter</td></tr><tr><td>L</td><td>0.005</td><td>0.005</td><td>length of the domain region</td></tr><tr><td>X<sub>I</sub></td><td>0</td><td></td><td>position of the interface</td></tr></tbody></table>

<table><caption>Table 3: Case 1: Computational and schemes parameters for simulation</caption><thead><tr><th>N<sub>p</sub></th><th>Δt</th></tr></thead><tbody><tr><td>2 × 10<sup>6</sup></td><td>0.01</td></tr><tr><td>θ<sub>g</sub></td><td>A</td><td>B</td><td>Ã</td><td>B̃</td></tr><tr><td>-0.4286</td><td>0.4286</td><td>0.5714</td><td>-0.4286</td><td>1.4286</td></tr></tbody></table>

We test the algorithms' performance by comparing the analytical ($q_{Lref}$) in (7)-(17) and numerical ($q_{Lnum}$) density transition function. As such, we verify whether the numerical schemes
accurately predict the particles distribution across the computational domain. Our estimator for the quality is the mean square error

$$RMSE = \sqrt{\frac{\sum_{i=1}^{N_a} (q_{Lref}(t, a_i) - q_{Lnum}(t, a_i))^2}{N_a}} \quad (46)$$

The comparison between the analytical and the numerical $q_L$, the latter obtained from the GUM, GSBM, and GHYMLA algorithms, is displayed in Figure 4 for the thermal point source initially located in $R_1$ and $R_2$.

## 4.2 Results

GUM and GSBM algorithms provide a good match between $q_{Lnum}$ and $q_{Lref}$ while GHYMLA algorithm shows a substantial discrepancy, especially for the case with the positive source (see Figure 4 (f)), being the thermal source closer to the discontinuity interface. Compared to the other two algorithms, GHYMLA requires a smaller time step to approximate the particles’ distribution with about the same accuracy due to the procedure used to estimate the first hitting time (see Appendix C.2). The RMSE between the analytical and the numerical $q_L$ is computed for both the initial positions of the point source (i.e., positive or negative source). The RMSE values at time $t = 0.5$ obtained from the GUM, GSBM and GHYMLA algorithms are reported in Table 4 along with the average computational time $t_{cm}$. GHYMLA is associated with RMSE values that are larger if compared with those obtained with GSBM and GUM, for the considered time level. Notably such performance difference increases in the case where the initial source is close to the interface (compare results between the case $X_0 = 2 \times 10^{-5}$ and $X_0 = -5 \times 10^{-5}$).

The average computational time differs significantly when different methods are compared. Due to the procedure used to calculate the first hitting time (see Appendix C.1), the GSBM algorithm is the most computationally expensive, GHYMLA and GUM being associated with similar computational efforts.
<figure><img src="image_6.png" /><figcaption>Figure 4: Case 1: Comparison between the analytical (black dashed line) and the numerical (red line) $q_L$ for a thermal point source initially located at $X_0 = -5 \times 10^{-5}$ (a,c,e) and $X_0 = 2 \times 10^{-5}$ (b,d,f) for the GUM (first row), GSBM (second row) and GHYMLA (third row) algorithm in the time range between 0.1 and 0.5 s with an interval of 0.1.</figcaption></figure>
<table><caption>Table 4: RMSE between the analytical and the numerical $q_L$ at time $t = 0.5$ and the average computational time $t_{cm}$, for a point source initially located at $X_0 = -5 \times 10^{-5}$ and $X_0 = 2 \times 10^{-5}$ for the GUM, GSBM and GHYMLA algorithms.</caption><thead><tr><th rowspan="2"></th><th colspan="2">$X_0 = -5 \times 10^{-5}$</th><th colspan="2">$X_0 = 2 \times 10^{-5}$</th></tr><tr><th>RMSE</th><th>$t_{cm}$</th><th>RMSE</th><th>$t_{cm}$</th></tr></thead><tbody><tr><td>GUM</td><td>12.122</td><td>0.1292</td><td>11.982</td><td>0.1379</td></tr><tr><td>GSBM</td><td>12.462</td><td>133.57</td><td>12.053</td><td>141.93</td></tr><tr><td>GHYMLA</td><td>14.401</td><td>0.0986</td><td>25.452</td><td>0.1032</td></tr></tbody></table>

We test the performance of the algorithms by varying the time step. RMSE between $q_{Lnum}$ and $q_{Lref}$ is calculated for $\Delta t$ spanning from 0.01 to 0.5 and shown in Figure 5. Over the analyzed range of $\Delta t$, the GSBM exhibits the lowest RMSE values and displays negligible sensitivity of the error with respect to $\Delta t$. The GUM algorithm is close to the result attained by GSBM for $\Delta t \le 0.05$ and shows a nonlinear increase (in log-log) for larger time steps. A direct linear convergence of the order of $\approx \frac{1}{2}$ is found for the GHYMLA. GHYMLA is the least accurate algorithm for the time step values analyzed.

Based on the RMSE and the computational time $t_{cm}$, the GUM algorithm offers the best balance between accuracy and computational load.

<figure><img src="image_7.png" /><figcaption>Figure 5: Case 1: RMSE between the analytical $q_{Lref}$ and the numerical $q_{Lnum}$ distribution obtained from the GUM (red circles), GSBM (blue squares) and GHYMLA (black diamonds) for different $\Delta t$ at time $t = 0.5$ and the thermal source at $X_0 = -5 \times 10^{-5}$.</figcaption></figure>
# 5 Conclusions

This paper provides a methodological extension that can be employed to approximate transport problems in heterogeneous domains adopting a Lagrangian approach.

Upon relying on the method of images, we analytically derived a closed-form equation for the fundamental solution $q(t, x, y)$ for diffusive transport in a discontinuous medium with generalized interface conditions at the discontinuity. We prove that assuming mass conservation, $q(t, x, \cdot)$ may be interpreted as a density of the position of a particle initially at $x$ at time $t$, leading to random walk methods. Based on these analytical expressions we propose a natural generalization of the Uffink, HYMLA and SBM methods, here referred as GUM, GHYMLA and GSBM. The initial problem (1)-(4) depends on four parameters $\Theta := (D_1, D_2, \nu_1/\nu_2, \kappa_1/\kappa_2) \in \mathbb{R}_+^4$. Mass conservation implies that the GUM, GSBM and GHYMLA methods can be applied when $\Theta$ belongs to a 3-dimensional manifold of $\mathbb{R}_+^4$, as it imposes one constraint (the original Uffink and HYMLA methods can be applied when $\Theta$ belongs to a 2-dimensional manifold as $\nu_i, \kappa_i$ are determined by the pair $(D_1, D_2)$). The presented algorithms are used to simulate diffusive transport in a bimaterial medium and benchmarked against an analytical solution.

The GUM algorithm offers the best compromise between accuracy and computational cost. Under the same time step, the GHYMLA algorithm exhibits the lowest accuracy in approximating the exact distribution of particles over the domain and a comparable computational time as GUM. Although exact, the GSBM method requires the simulation of complex distributions, which hinders its computational time. A significant impact of the time step on the GUM ($\Delta t > 0.05$) and GHYMLA is observed, while the GSBM is not affected by $\Delta t$.

Our methods can be applied for any value assumed by the coefficients $\nu_i$ and $\kappa_i$ describing the jump of the solution and its gradient at the interface under a condition of mass conservation linking these parameters. This covers a wide range of differential problems representing diverse physical transport phenomena (see Appendix D). An extension of the methods proposed here to remove the constraint of the mass conservation will be presented in [2]. In particular, we give a practical application to the Fourier law.

# References

[1] T. Appuhamillage, V. Bokil, E. Thomann, E. Waymire, and B. Wood. Solute transport across an interface: A fickian theory for skewness in breakthrough curves. *Water Resources Research*, 46(7), 2010.
[2] E. Baioni, A. Lejay, G. Pichot, and G. M. Porta. Random walk modeling of conductive heat transport in discontinuous media. working paper or preprint, July 2023.

[3] E. Baioni, M. Mousavi Nezhad, G. M. Porta, and A. Guadagnini. Modeling solute transport and mixing in heterogeneous porous media under turbulent flow conditions. *Physics of Fluids*, 33(10), 2021.

[4] M. Bechtold, J. Vanderborght, O. Ippisch, and H. Vereecken. Efficient random walk particle tracking algorithm for advective-dispersive transport in media with discontinuous dispersion coefficients and water contents. *Water Resources Research*, 47(10), 2011.

[5] D. A. Benson, T. Aquino, D. Bolster, N. Engdahl, C. V. Henri, and D. Fernandez-Garcia. A comparison of eulerian and lagrangian transport and non-linear reaction algorithms. *Advances in Water Resources*, 99:15–37, 2017.

[6] B. Berkowitz, A. Cortis, I. Dror, and H. Scher. Laboratory experiments on dispersive transport across interfaces: The role of flow direction. *Water Resources Research*, 45(2), 2009.

[7] T. G. Bhaskar and R. Kumar. Analyticity of semigroup generated by a class of differential operators with interface. *Nonlinear Analysis*, 39(6, Ser. A: Theory Methods):779–791, 2000.

[8] S. Chapman. On the brownian displacements and thermal diffusion of grains suspended in a non-uniform fluid. *Proceedings of the Royal Society of London. Series A, Containing Papers of a Mathematical and Physical Character*, 119(781):34–54, May 1928.

[9] C. Cordes, H. Daniels, and G. Rouvé. New very efficient algorithm for particle tracking in layered aquifers. In *The 2nd International Conference on Computer Methods in Water Resources, Marrakesh, Morocco, 02/20-22/91*, pages 41–55, 1991.

[10] W. Feller. *An introduction to probability theory and its applications, Vol II*. John Wiley & Sons, Inc., New York-London-Sydney, second edition, 1971.

[11] C. V. Henri and E. Diamantopoulos. Unsaturated transport modeling: Random-walk particle-tracking as a numerical-dispersion free and efficient alternative to eulerian methods. *Journal of Advances in Modeling Earth Systems*, 14(9):e2021MS002812, 2022.

[12] H. Hoteit, R. Mose, A. Younes, F. Lehmann, and P. Ackerer. Three-dimensional modeling of mass transfer in porous media using the mixed hybrid finite elements and the random-walk methods. *Mathematical Geology*, 34(4):435–456, 2002.
[13] E. M. LaBolle, G. E. Fogg, and A. F. Tompson. Random-walk simulation of transport in heterogeneous porous media: Local mass-conservation problem and implementation methods. *Water Resources Research*, 32(3):583–593, 1996.

[14] E. M. Labolle, J. Quastel, G. E. Fogg, and J. Gravner. Diffusion processes in composite porous media and their numerical integration by random walks: Generalized stochastic differential equations with discontinuous coefficients. *Water Resources Research*, 36:651–662, 2000.

[15] A. Lejay. On the constructions of the skew Brownian motion. *Probability Surveys*, 3:413–466, 2006.

[16] A. Lejay and G. Pichot. Simulating diffusion processes in discontinuous media: a numerical scheme with constant time steps. *Journal of Computational Physics*, 231(21):7299–7314, 2012.

[17] A. Lejay and G. Pichot. Simulating diffusion processes in discontinuous media: Benchmark tests. *Journal of Computational Physics*, 314:384–413, 2016.

[18] D.-H. Lim. Numerical study of nuclide migration in a nonuniform horizontal flow field of a high-level radioactive waste repository with multiple canisters. *Nuclear Technology*, 156(2):222–245, 2006.

[19] Y. Maruyama. Random walk to describe diffusion phenomena in three-dimensional discontinuous media: Step-balance and fictitious-velocity corrections. *Physical Review E*, 96(3):032135, 2017.

[20] Y. Maruyama. Random walk to describe diffusion phenomena in three-dimensional discontinuous media: Extension to anisotropic phases and nonplane interfaces. *Physical Review E*, 101(2):023307, 2020.

[21] H. Oukili, R. Ababou, G. Debenest, and B. Noetinger. Random walks with negative particles for discontinuous diffusion and porosity. *Journal of Computational Physics*, 396:687–701, 2019.

[22] J. M. Ramirez, E. A. Thomann, E. C. Waymire, J. Chastanet, and B. D. Wood. A note on the theoretical foundations of particle tracking methods in heterogeneous porous media. *Water Resources Research*, 44(1), 2008.

[23] P. Salamon, D. Fernandez-Garcia, and J. Gómez-Hernández. Modeling tracer transport at the made site: The importance of heterogeneity. *Water Resources Research*, 43(8), 2007.
[24] M. J. Schmidt, N. B. Engdahl, S. D. Pankavich, and D. Bolster. A mass-transfer particle-tracking method for simulating transport with discontinuous diffusion coefficients. *Advances in Water Resources*, 140:103577, 2020.

[25] K. Semra, P. Ackerer, and R. Mosé. Three dimensional groundwater quality modelling in heterogeneous media. *WIT Transactions on Ecology and the Environment*, 2, 1970.

[26] G. Uffink. A random walk method for the simulation of macrodispersion in a stratified aquifer. *Relation of Groundwater Quantity and Quality*, pages 103–114, 1983.

[27] G. Uffink. Analysis of dispersion by the random walk method. *Ph.D. Dissertation*, Delft University of Technology, 1990.

[28] J. Walsh. A diffusion with discontinuous local time. In *Temps locaux*, volume 52–53 of *Astérisques*, pages 37–45. Société Mathématique de France, 1978.

[29] M. T. Wereide. La diffusion d'une solution dont la concentration et la température sont variables. *Annales de Physique*, 9(2):67–83, 1914.

[30] J. A. Yan. A formula for densities of transition functions. In *Séminaire de Probabilités*, XXII, volume 1321 of *Lecture Notes in Mathematics*, pages 92–100. Springer, Berlin, 1988.

[31] Y. Zhang, E. M. LaBolle, and K. Pohlmann. Monte carlo simulation of superdiffusion and subdiffusion in macroscopically heterogeneous media. *Water Resources Research*, 45(10), 2009.

# Appendices

## A A convenient transform

Let $f$ be a function with a condition of type

$$ \kappa_1 f'(0^-) = \kappa_2 f'(0^+) $$

for some $\kappa_1, \kappa_2 > 0$. Such a condition depends actually on the ratio $\kappa_1/\kappa_2$. We introduce the Möbius transform

$$ \begin{aligned} \mathbf{M}: (0, +\infty) &\to (-1, 1) \\ z &\mapsto \frac{1-z}{1+z}. \end{aligned} $$
This function has the remarkable property that $M^{-1} = M$. When

$$ \gamma := M \left( \frac{\kappa_1}{\kappa_2} \right) = \frac{\kappa_2 - \kappa_1}{\kappa_2 + \kappa_1}, $$

we obtain that

$$ \frac{\kappa_1}{\kappa_2} = \frac{1 - \gamma}{1 + \gamma} $$

so that (A) becomes

$$ (1 - \gamma) f'(0^{-}) = (1 + \gamma) f'(0^{+}). \qquad (47) $$

# B Mass conservation

We check that the mass conservation property given in (24) holds under the necessary and sufficient conditions that $A + B = \tilde{A} + \tilde{B} = 1$, where the constants $A, B, \tilde{A}$ and $\tilde{B}$ defined in Sections 3.1.1 and 3.1.2.

## B.1 Negative source $x < 0$

Integrating the densities $p_1$ and $p_2$,

$$ \int_{0}^{+\infty} p_2(t, \beta x, y) \, dy = \int_{0}^{+\infty} \beta p_2(t, \beta x, \beta y) \, dy = \int_{0}^{+\infty} p_1(t, x, y) \, dy. \qquad (48) $$

Therefore,

$$ \int_{-\infty}^{+\infty} q(t, x, y) \, dy = \int_{-\infty}^{0} p_1(t, x, y) \, dy + A \int_{-\infty}^{0} p_1(t, -x, y) \, dy + B \int_{0}^{+\infty} p_1(t, x, y) \, dy. \qquad (49) $$

As $p_1(t, -x, y) = p_1(t, x, -y)$, it follows that

$$ \int_{-\infty}^{0} p_1(t, -x, y) \, dy = \int_{0}^{+\infty} p_1(t, x, y) \, dy. \qquad (50) $$

It follows that

$$ \int_{-\infty}^{+\infty} q(t, x, y) \, dy = \int_{-\infty}^{0} p_1(t, x, y) \, dy + (A + B) \int_{0}^{+\infty} p_1(t, x, y) \, dy. \qquad (51) $$

Hence, the mass is conserved if and only if $A + B = 1$.
## B.2 Positive source $x > 0$

Similarly to the case $x < 0$ we consider

$$ \int_{-\infty}^{0} p_{1}(t, \tilde{\beta}x, y) dy = \int_{-\infty}^{0} \tilde{\beta} p_{1}(t, \tilde{\beta}x, \tilde{\beta}y) dy = \int_{-\infty}^{0} p_{2}(t, x, y) dy. \quad (52) $$

It follows that

$$ \int_{-\infty}^{+\infty} q(t, x, y) dy = \tilde{B} \int_{-\infty}^{0} p_{2}(t, x, y) dy + A \int_{0}^{+\infty} p_{2}(t, x, y) dy + \tilde{A} \int_{0}^{+\infty} p_{2}(t, -x, y) dy. \quad (53) $$

As $p_2(t, -x, y) = p_2(t, x, -y)$, we get

$$ \int_{0}^{+\infty} p_{1}(t, -x, y) dy = \int_{-\infty}^{0} p_{2}(t, x, y) dy. \quad (54) $$

Therefore,

$$ \int_{-\infty}^{+\infty} q(t, x, y) dy = (\tilde{A} + \tilde{B}) \int_{-\infty}^{0} p_{2}(t, x, y) dy + \int_{0}^{+\infty} p_{2}(t, x, y) dy. \quad (55) $$

Hence, the mass is conserved if and only if $\tilde{A} + \tilde{B} = 1$.

# C GSBM and GHYMLA algorithms

## C.1 GSBM algorithm

Similarly to the SBM algorithm [16, 17], the GSBM consists of two main steps to advance the particle's position in a time step $\Delta t$: (Step 1) computation of the first hitting time $\tau$ according to the *ExactHittingTime* algorithm [17], which is grounded on the Brownian Bridge concept, and (Step 2) the computation of the final particle's position based on (Step 1). If $\tau > \Delta t$ the particle initially at $x$ does not reach the discontinuity interface during the time step $\Delta t$ and the displacement is computed as $\epsilon\sqrt{2D(x)\Delta t}$ where $\epsilon$ is a random variate from a Gaussian distribution $N(0, 1)$. When $\tau < \Delta t$ the particle's location at time $t + \Delta t$ corresponds to the distance the particle travels from $X_I = 0$ during the interval $\Delta t - \tau$.

The GSBM algorithm is reported in Algorithm 2.
**Algorithm 2: GSBM algorithm**

**Input:** An initial position $X_t = x$ and a time step $\Delta t > 0$ with an interface at 0.

**Output:** Particle's position $X_{t+\Delta t}$ according to GSBM algorithm.

/* Compute the first hitting time $\tau$ */

Generate a random variate $\xi \sim N(0, 1)$

Set $z = x/\sqrt{2D(x)}$

Set $y = z + \xi\sqrt{\Delta t}$

if $zy<0$ then

/* The interface is crossed */

Generate a random variate $\xi \sim IG(|z|/|y|, z^2/\Delta t)$

Set $\tau = \Delta t \xi/(\xi+1)$

else

/* Check if the interface has been crossed */

Generate a random variate $U \sim U(0, 1)$

if $U < \exp(-2zy/\Delta t)$ then

/* The interface has been crossed */

Generate a random variate $\xi \sim IG(|z|/|y|, z^2/\Delta t)$

Set $\tau = \Delta t \xi/(\xi+1)$

else

/* The interface has not been crossed */

Set $\tau = t + \Delta t$

end

end

/* Compute the final particle's position */

if $\tau < \Delta t$ then

/* The interface has been crossed */

Generate a random variate $U \sim U(0, 1)$

Generate a random variate $\epsilon \sim N(0, 1)$

if $U < \frac{1+\theta_q}{2}$ then

return $\sqrt{2D_2(\Delta t - \tau)}|\epsilon|$

else

return $-\sqrt{2D_1(\Delta t - \tau)}|\epsilon|$

end

else

/* The interface has not been crossed */

return $y\sqrt{2D(y)}$

end
## C.2 GHYMLA algorithm

The GHYMLA algorithm presented in this paper and reported in Algorithm 3 is a generalized version of the original HYMLA algorithm proposed by Hoteit, Younes, Mose, Lehmann and Ackerer [12].

GHYMLA is “a splitting up jump” algorithm comprising two steps: (Step 1) estimating the time the particle needs to reach the discontinuity (first hitting time $\tau$), and (Step 2) determining the final particle’s position: if the particle arrives at $X_I$ during the time step $\Delta t$, i.e., $\tau < \Delta t$, the particle’s location at time $t + \Delta t$ corresponds to displacement during the time $\Delta t_{rem} = \Delta t - \tau$ starting from the $X_I$, if $\tau > \Delta t$ the particle moves according to a standard random walk.

The first hitting time is estimated considering a linear relation between the position of the particle at time $t$ ($x$) and $t + \Delta t$ ($y$), the latter being quantified as $(2U - 1)\sqrt{6D\Delta t}$ where $U$ is a random variate from a uniform distribution $U(0, 1)$. According to the LinearHittingTime algorithm [17], $\tau$ results

$$ \tau = \frac{|x|}{|x| + |y|} \Delta t. \qquad (56) $$

Once the particle reaches the discontinuity, the probability to move towards the region $y < 0$ ($R_1$) or $y > 0$ ($R_2$) starting from 0 (see Figure 6) is quantified as [12]

$$ \mathbb{P}_1 = \frac{1 - \theta_g}{2} \quad \text{and} \quad \mathbb{P}_2 = 1 - \mathbb{P}_1 = \frac{1 + \theta_g}{2}. \qquad (57) $$

This consists in replacing the properly renormalized half-Gaussian random variates when the particle is at 0 in the GSBM scheme by uniform random variates. The values of the pair $(\mathbb{P}_1, \mathbb{P}_2)$ in (57) corresponds to the one of the pair $(\mathbb{P}_-, \mathbb{P}_+)$ in (30).

<figure><img src="image_8.png" /><figcaption>Figure 6: Density for a particle to move to the region $y < 0$ ($R_1$) or $y > 0$ ($R_2$) starting from 0.</figcaption></figure>

We recover the HYMLA algorithm when the density of the particle from 0 is uniformly distributed over $[-\sqrt{6D_1\Delta t_{rem}}, \sqrt{6D_2\Delta t_{rem}}]$, that is when

$$ \frac{\mathbb{P}_1}{\sqrt{6D_1\Delta t_{rem}}} = \frac{\mathbb{P}_2}{\sqrt{6D_2\Delta t_{rem}}}. $$
**Algorithm 3: GHYMLA algorithm**

**Input:** An initial position $X_t = x$ and a time step $\Delta t > 0$ with an interface at $x_I = 0$.

**Output:** Particle's position $X_{t+\Delta t}$ at time $t + \Delta t$ according to GHYMLA algorithm.

Compute $H_1 = \sqrt{6D_1\Delta t}$ and $H_2 = \sqrt{6D_2\Delta t}$

Generate a random variate $V \sim U(0, 1)$

if $x < 0$ then

$$y = x + (2V - 1)H_1$$

else

$$y = x + (2V - 1)H_2$$

end

/* Compute the first hitting time $\tau$ */

if $zy < 0$ then

$$ \left| \begin{array}{l} \text{/* The interface has been crossed} \\ \tau = \frac{|x|}{|x|+|y|}\Delta t \end{array} \right. \quad \text{*/} $$

else

$$ \left| \begin{array}{l} \text{/* The interface has not been crossed} \\ \tau = \Delta t \end{array} \right. \quad \text{*/} $$

end

/* Compute the particle's final position */

if $\tau < \Delta t$ then

$$ \left| \begin{array}{l} \text{/* The interface is crossed: biased step} \\ \text{/* Time splitting} \end{array} \right. \quad \text{*/} $$

Compute $\Delta t_{\text{rem}} = \Delta t - \tau$

Generate a random variate $U \sim U(0, 1)$

if $U < \frac{1-\theta_q}{2}$ then

$$ \left| \begin{array}{l} \text{Compute } x_L = -\sqrt{6D_1\Delta t_{\text{rem}}} \\ \text{Set } x_R = 0 \end{array} \right. $$

else

$$ \left| \begin{array}{l} \text{Set } x_L = 0 \\ \text{Compute } x_R = \sqrt{6D_2\Delta t_{\text{rem}}} \end{array} \right. $$

end

Generate a random variate $U \sim U(0, 1)$

return $x_L + (x_R - x_L)U$

else

$$ \left| \begin{array}{l} \text{/* The interface is not crossed: uniform step} \\ \text{return } y \end{array} \right. \quad \text{*/} $$

end
Since $P_1 + P_2 = 1$, this is equivalent to

$$ \theta_g = \frac{\sqrt{D_2} - \sqrt{D_1}}{\sqrt{D_1} + \sqrt{D_2}}. $$

This is the same condition as the one to apply the original Uffink method and correspond to the Fick law (see Section D below).

# D Examples for physically relevant equations

(I) The Fick or Darcy law, corresponds to the divergence-form operator $\nabla_y(D(y)\nabla_y\cdot)$, so that $y \mapsto q(t, x, y)$ and $y \mapsto D(y)\nabla_y q(t, x, y)$ shall be continuous. Therefore

$$ \nu_1 = \nu_2 = 1 \text{ and } \kappa_1 = D_1, \kappa_2 = D_2. $$

Then the mass is conserved, and

$$ \Gamma = \frac{\sqrt{D_1}}{\sqrt{D_2}}, \quad P_+ = \frac{1 + \theta_g}{2} = \frac{\sqrt{D_2}}{\sqrt{D_2} + \sqrt{D_1}}. $$

(II) The Chapman equation [8] corresponds to the operator $\Delta_y(D(y)\cdot)$, meaning that

$$ \nu_1 = D_1, \nu_2 = D_2, \kappa_1 = D_1, \text{ and } \kappa_2 = D_2. $$

since $y \mapsto D(y)q(t, x, y)$ and its gradient shall be continuous. Hence, the mass is conserved, and

$$ \Gamma = \frac{\sqrt{D_2}}{\sqrt{D_1}} \text{ and } P_+ = \frac{1 + \theta_g}{2} = \frac{\sqrt{D_2}}{\sqrt{D_2} + \sqrt{D_1}}. $$

The stochastic process associated with this operator is the Itô stochastic differential equation.

(III) The Wereide equation [29] corresponds to the operator $\nabla_y(\sqrt{D(y)}\nabla_y(\sqrt{D(y)}\cdot)$ in which the diffusion is due to osmotic pressure and to particles moving irregularly in the media. In this case,

$$ \nu_1 = \sqrt{D_1}, \nu_2 = \sqrt{D_2}, \kappa_1 = D_1, \kappa_2 = D_2. $$

The mass is conserved, and

$$ \Gamma = \frac{D_2}{D_1} \text{ and } P_+ = \frac{1 + \gamma_g}{2} = \frac{D_1}{D_1 + D_2}. $$

Only the first case may be treated with the HYMLA and Uffink methods. The two other cases require their generalized versions.