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

HAL
open science

# Random walk modeling of conductive heat transport in discontinuous media

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

► To cite this version:

Elisa Baioni, Antoine Lejay, Géraldine Pichot, Giovanni Michele Porta. Random walk modeling of conductive heat transport in discontinuous media. 2023. [https://arxiv.org/abs/2308.04652](https://arxiv.org/abs/2308.04652)

HAL Id: hal-04166562

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

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
Random walk modeling
of conductive heat transport
in discontinuous media

Elisa Baioni,\* Antoine Lejay,† Géraldine Pichot,‡ Giovanni Michele Porta§

September 25, 2024

Published in *Transport in Porous Media*, DOI: 10.1007/s11242-024-02132-6

Abstract

We consider heat transport within a discontinuous domain by relying on the modeling approach proposed by Baioni et al. Such approach has been specifically designed to address diffusive processes in media with discontinuous physical properties and generalized boundary conditions at the discontinuities. Three algorithms are here applied to estimate the conductive heat transport in a bimaterial medium. The algorithms undergo testing using two test cases that share the same computational domain but differ in terms of their initial conditions. According to the numerical results all the algorithms ensure the conservation of thermal energy and preserve thermal equilibrium under steady state conditions. The Generalized Uffink Method (GUM) and Generalized HYMLA (GHYMLA) demonstrate sensitivity to the choice of the time step, whereas the Generalized Skew Brownian Motion (GSBM) appears to be unaffected by the value of $\Delta t$. The GUM algorithm presents an optimal trade-off between accuracy and computational time.

**Keywords:** Discountinuous media, thermal diffusivity, conductive heat transport, numerical simulation, random walk methods

# 1 Introduction

Discontinuous media, characterized by irregular interfaces and varying material properties, are encountered in a wide range of systems, including porous materials, composite structures, and multiphase systems. Heat transfer across different media is a phenomenon that affects numerous environmental [24, 28] and engineering processes [6, 9, 14, 22, 32] and stands as a vital component of clean and sustainable energy systems. For example, porous media play a crucial role in heat exchange within electronic cooling, heat exchangers, and various thermal systems. As discussed

<footnote>\*Department of Civil and Environmental Engineering, Politecnico di Milano, Piazza Leonardo da Vinci, Milan, 20133, Italy. elisa.baioni@polimi.it</footnote>

<footnote>†CNRS, Inria, IECL, Université de Lorraine, Boulevard des Aiguillettes, Nancy, F-54000, France. Antoine.Lejay@univ-lorraine.fr</footnote>

<footnote>‡Inria, 2 rue Simone Iff, Paris, 75589, France & CERMICS (ENPC), Université Paris-Est, 6 et 8 av. Blaise Pascal, Marne-la-Vallée Cedex 2, 77455, France, geraldine.pichot@inria.fr</footnote>

<footnote>§Department of Civil and Environmental Engineering, Politecnico di Milano, Piazza Leonardo da Vinci, Milan, 20133, Italy. giovanni.porta@polimi.it</footnote>
Random walk modeling of conductive heat transport
in discontinuous media

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

¹Politecnico di Milano, Department of Civil and Environmental Engineering,
20133 Milan (IT); ORCID:0000-0001-8146-7435

²Université de Lorraine, CNRS, Inria, IECL, F-54000 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

⁴Politecnico di Milano, Department of Civil and Environmental Engineering,
20133 Milan (IT); ORCID:0000-0002-0636-373X

Abstract

We consider heat transport within a discontinuous domain by relying on the modeling approach proposed by Baioni et al. Such approach has been specifically designed to address diffusive processes in media with discontinuous physical properties and generalized boundary conditions at the discontinuities. Three algorithms are here applied to estimate the conductive heat transport in a bimaterial medium. The algorithms undergo testing using two test cases that share the same computational domain but differ in terms of their initial conditions. According to the numerical results all the algorithms ensure the conservation of thermal energy and preserve thermal equilibrium under steady state conditions. The Generalized Uffink Method (GUM) and Generalized HYMLA (GHYMLA) demonstrate sensitivity to the choice of the time step, whereas the Generalized Skew Brownian Motion (GSBM) appears to be unaffected by the value of $\Delta t$. The GUM algorithm presents an optimal trade-off between accuracy and computational time.

**Keywords–** Discountinuous media, thermal diffusivity, conductive heat transport, numerical simulation, random walk methods
in various studies [13,26], the existence of abrupt changes in material properties and/or physical conditions within a medium can introduce significant challenges when modeling heat transfer phenomena.

Conductive heat transfer in layered composite materials and functionally graded materials [1] has been widely addressed in the literature, where analytical solutions have been presented [2,7] for simple geometries and a prescribed set of initial and boundary conditions. To overcome these limiting assumptions numerous literature works deal with numerical simulation tools.

Porous media are inherently characterized by the presence of material discontinuities, typically at the interface between the solid matrix and the pore space. Diverse approximation methods have been proposed in the literature to deal with transport in such conditions. Direct solvers feature the application of computational fluid dynamics tools [22], these methods being limited to small size samples due to computational constraints. The simulation of heat transfer has also been addressed through dedicated pore network models (PNM). The pore network model involves the conceptualization of the geometry of the porous medium through an interconnected network of simplified pores and throats [31]. This approach effectively reduces the computational demand compared to direct numerical simulation. However, the accuracy of the model is constrained by its ability to capture the intricate topological details of the original pore structure. Consequently, the PNM may struggle to accurately model scenarios involving heterogeneous temperature distributions. Koch et al. [15–17] introduced the dual network model (DNM), which simultaneously simulates coupled single-phase flow and energy transport in porous media, enabling the accurate representation of heterogeneous temperature distributions in both the void space and solid matrix. Stochastic methods have also been developed starting from a general Monte Carlo framework [8, 11, 27]. Stochastic Monte Carlo approaches have been proposed to simulate heat transport in complex media by obtaining temperature through a random sampling of pathways, typically resulting from a discretized version of Fourier’s law [8,27]. These methods can become computationally demanding, particularly for multi-dimensional problems.

This study addresses heat transport in a one-dimensional discontinuous medium relying on the random walk approaches introduced by [4]. These latter enable the resolution of diffusive processes in domains characterized by discontinuous physical properties and generalized interface conditions. Random walk solvers are a particular class of Monte Carlo method, that have been typically applied to resolve advective-diffusive solute transport. Appropriate stochastic processes may be applied in this context to address various type of diffusive phenomena [18, 19, 23, 29].

Conductive heat transfer typically follows a diffusive dynamic and has been resolved using particle tracking approaches in homogeneous media [11]. Random walk particle tracking approaches have proven useful as simulation methods and have facilitated the development of conceptual and mathematical models to upscale pore-scale properties [25,29]. For diffusive mass transfer several algorithms can be applied to deal with the presence of spatially discontinuous diffusion coefficients [20]. However, to the best of our knowledge, random walk methods directly based on a diffusive stochastic process in discontinuous media have not been formally proposed for heat transfer in existing literature works. Drawing an analogy with the case of transport of dissolved mass, heat transfer can be resolved assuming energy conservation. However, Fourier’s law implies a diffusive process where temperature is the variable of interest, which is not conserved. Notably, the coefficient governing heat transfer by conduction within a material (often called thermal diffusivity) is different than the conductivity coefficient, which is conversely applied to guarantee heat flux continuity across a permeable interface. This prevents direct applications of existing random walk algorithms to heat transfer problems in discontinuous media.

In this work we address the formal derivation of random walk algorithms for heat transfer in discontinuous media. Here the differential problem (Fourier's law) and the related interface conditions at material discontinuities are interpreted in terms of the underlying stochastic
# 1 Introduction

Discontinuous media, characterized by irregular interfaces and varying material properties, are encountered in a wide range of systems, including porous materials, composite structures, and multiphase systems. Heat transfer across different media is a phenomenon that affects numerous environmental [19, 24] and engineering processes [12, 7, 4, 18, 28] and stands as a vital component of clean and sustainable energy systems. For example, porous media play a crucial role in heat exchange within electronic cooling, heat exchangers, and various thermal systems. As discussed in various studies [11, 22, 21], the existence of abrupt changes in material properties and/or physical conditions within a medium can introduce significant challenges when modeling heat transfer phenomena.

Conductive heat transfer in layered composite materials has been widely addressed in the literature, where analytical solutions have been presented [5, 1] for simple geometries and a prescribed set of initial and boundary conditions. To overcome these limiting assumptions numerous literature works deal with numerical simulation tools. Numerical analysis of heat transfer processes can be performed by direct numerical simulation based on computational fluid dynamics tools (e.g., [18]), these methods being limited to small size samples due to computational constraints. Alternative approaches are based on probabilistic simulations based on Monte Carlo approaches [9, 6, 27, 23], or on pore network modeling [14, 13, 15], that allow analyzing larger systems by conceptualizing the geometry through a interconnected network of simplified pores and throats.

This study addresses heat transport in a discontinuous medium relying on the random walk approaches introduced in [2]. These latter enable the resolution of diffusive processes in domains characterized by discontinuous physical properties and generalized interface conditions. Random walk particle tracking approaches have been widely employed to simulate advective-diffusive mass transfer in porous domains. In this context these class of Lagrangian methods have proven useful as simulation methods and have facilitated the development of conceptual and mathematical models to upscale pore-scale properties [20, 25]. For diffusive mass transfer several algorithms can be applied to deal with the presence of spatially discontinuous diffusion coefficients [17]. Yet, the same methodologies cannot be applied to heat transfer, because of the inherent difference in the problem parameterization between heat and mass transfer. Notably, the coefficient driving heat transfer by conduction within a material (often called thermal diffusivity) is different than the conductivity coefficient, which is conversely applied to guarantee heat flux continuity across an interface. The applicability of classical random walk methods to heat transfer processes is thus strictly limited to homogeneous media. A general theory and solution algorithms are developed in [2] to deal with generalized interface conditions across an interface. In this paper we explicitly show how such theory can be leveraged to
processes. For this reason, our random walk method differs from others Monte Carlo approaches that have been defined starting from different approximations of Fourier's law, i.e. typically employing MC methods as a stochastic solver for a discretized version of the partial differential equation [8,27]. To provide a numerical approximation, we employ a set of three random walk algorithms recently introduced in [4] for a wide class of diffusive processes. These methodologies are here specifically adapted to and benchmarked for heat transfer.

A general theory and solution algorithms are developed by [4] to deal with generalized permeable interface conditions across an interface. In this paper we explicitly show how such theory can be leveraged to simulate conductive heat transport across materials characterized by discontinuous properties. To overcome the above mentioned limitations of existing methods in approximating Fourier's law, which governs heat transport, we have proposed a new formulation through a convenient transformation implementing a Lagrangian method with a weight associated to each particle. Ultimately the developed approaches provide a critical extension of available Lagrangian methods and pave the way towards the application of such methods to a new class of problems. We employ two test cases dealing with heat transfer in a bimaterial domain and benchmarked against an analytical solution. Two test cases, each involving the same computational domain but different initial conditions, are presented to evaluate the performance of the algorithms in estimating heat transport in discontinuous media.

**The paper is structured as:** Section 2 presents the application of the three generalized interface algorithms, i.e. the generalized Uffink, HYMLA and Skew Brownian Motion algorithms, here labeled as GUM, GHYMLA and GSBM respectively. The benchmark tests used to evaluate the algorithms' performance are described in Section 3, as well as a discussion of the numerical results. The conclusions drawn from our study are reported in Section 4.

## 2 Heat transport in a bimaterial medium

We consider solving the Fourier's law [10] with an interface at 0 and an initial condition $\psi$:

$$ \phi(y)\partial_t T(t, y; \psi) = \mathcal{F}_y T(t, y; \psi), \quad (1) $$

$$ T(t, 0^{-}; \psi) = T(t, 0^{+}; \psi), \quad (2) $$

$$ k_1 \nabla T(t, 0^{-}; \psi) = k_2 \nabla T(t, 0^{+}; \psi), \quad (3) $$

$$ T(0, y; \psi) = \psi(y), \quad (4) $$

where $T(t, y; \psi)$ is the temperature at time $t$ and point $y$,

$$ \mathcal{F}_y f(y) := \nabla_y(k(y)\nabla_y f(y)) \quad \text{and} \quad \phi(y) := \rho(y)c(y). $$

Here, $\rho$, $c$ and $k$ are functions constant on each side of the permeable interface located at 0. For such functions $f$, we write $f(y) = f_1$ if $y < 0$ and $f_2$ if $y > 0$. The constants are the thermal conductivity $k_i$ [W m$^{-1}$ K$^{-1}$], the specific heat $c_i$ [J kg$^{-1}$ K$^{-1}$], and the density $\rho_i$ [kg m$^{-3}$]. All these constants are assumed to be positive.

### 2.1 Properties of heat transport

First, we consider the problem in free space, and later in a finite domain $[-L, L]$ with reflected boundary conditions at each side. The computations are the same in such a case.
simulate conductive heat transport across materials characterized by discontinuous properties.
To overcome the above mentioned limitations of existing methods in approximating Fourier's
law, which governs heat transport, we have proposed a new formulation through a convenient
transformation implementing a Lagrangian method with a weight associated to each particle.
Ultimately the developed approaches provide a critical extension of available Lagrangian
methods and pave the way towards the application of such methods to a new class of problems.
We employ two test cases dealing with heat transfer in a bimaterial domain and benchmarked
against an analytical solution. Two test cases, each involving the same computational domain
but different initial conditions, are presented to evaluate the performance of the algorithms in
estimating heat transport in discontinuous media.

**The paper is structured as:** Section 2 presents the application of the three generalized interface algorithms, i.e. the generalized Uffink, HYMLA and Skew Brownian Motion algorithms, here labeled as GUM, GHYMLA and GSBM respectively. The benchmark tests used to evaluate the algorithms' performance are described in Section 3, as well as a discussion of the numerical results. The conclusions drawn from our study are reported in Section 4.

## 2 Heat transport in a bimaterial medium

We start from the approach implemented in [2] to solve diffusive transport in discontinuous
media under generalized conditions at the discontinuity interface. Here we consider one-
dimensional thermal diffusion process through a random walk. The GUM, GHYMLA and
GSBM algorithms reported in the Supplementary Material are used to assess conductive heat
transfer in a bimaterial medium, as in Figure 1, 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 a diffusivity equal to $D_1$
and $D_2$, respectively. The theoretical formulation of the physical problem and the derivation
of the model constants for the case study are illustrated in the following.

### 2.1 Properties of heat transport

The GUM, GSBM and GHYMLA algorithms discussed in [2], solve diffusive problems in
discontinuous media with generalized boundary conditions at the discontinuity interface.

We now present how to apply these methods to a heat transport problem. Actually, this kind
of problem cannot be treated directly by applying GUM, GSBM or GHYMLA so that we
introduce a convenient transform. More precisely, we consider solving the Fourier's law [8]
Integrating (1) on the domain and using the Green formula,

$$ \begin{aligned} \partial_t \int_{-\infty}^{+\infty} \phi(y) T \, dy &= \int_{-\infty}^{0} \nabla_y(k(y)\nabla_y T) \, dy + \int_{0}^{+\infty} \nabla_y(k(y)\nabla_y T) \, dy \\ &= k_1 \nabla T(t, 0^{-}; \psi) - k_2 \nabla T(t, 0^{+}; \psi) = 0. \quad (5) \end{aligned} $$

Hence, the quantity $\int_{-\infty}^{+\infty} \phi(y)T \, dy$ is conserved over time. On the other hand we observe that,

$$ \begin{aligned} \partial_t \int_{-\infty}^{+\infty} T \, dy &= \int_{-\infty}^{0} \frac{1}{\phi_1} \nabla_y(k_1 \nabla_y T) \, dy + \int_{0}^{+\infty} \frac{1}{\phi_2} \nabla_y(k_2 \nabla_y T) \, dy \\ &= \frac{k_1}{\phi_1} \nabla T(t, 0^{-}; \psi) - \frac{k_2}{\phi_2} \nabla T(t, 0^{+}; \psi). \quad (6) \end{aligned} $$

Therefore, the quantity $\int_{-\infty}^{+\infty} T \, dy$ is not conserved over time.

As a consequence, a convenient transform, explained in the next paragraph, is needed to fall into the class of diffusive problems in discontinuous media with generalized boundary conditions as considered in [4].

## 2.2 Lagrangian approach

The design of a diffusive Lagrangian scheme is achieved by defining a fundamental solution $q(t, x, y)$ [m⁻¹] such that $q(t, x, y)$ $dy$ is the probability that a particle starting from $x$ ends in a small volume $dy$ around $y$ at time $t$. As practical implication of (5) we design a numerical scheme relying upon a Lagrangian approach where the thermal energy of the system is discretized by a set of particles, each one carrying a specific quantity. We relate $q$ and $T$ through the convenient transform

$$ T = \int \psi(x)q(t, x, y) \frac{\phi(x)}{\phi(y)} \, dx. \quad (7) $$

Here, $q$ is the family of solutions to

$$ \partial_t q(t, x, y) = \mathcal{L}_y q(t, x, y), \quad (8) $$

$$ \frac{1}{\phi_1} q(t, x, 0^{-}) = \frac{1}{\phi_2} q(t, x, 0^{+}), \quad (9) $$

$$ \alpha_1 \nabla q(t, x, 0^{-}) = \alpha_2 \nabla q(t, x, 0^{+}), \quad (10) $$

$$ \lim_{t \to 0} q(t, x, y) = \delta_x(y), \quad (11) $$

where

$$ \mathcal{L}_y f(y) = \nabla_y(\alpha(y)\nabla_y f(y)) \quad \text{with} \quad \alpha_i := \frac{k_i}{\rho_i c_i} \, [\text{W m}^2 \text{ J}^{-1}], \quad i = 1, 2, $$

for any suitable function $f$. The quantity $\alpha$ is the thermal diffusivity.

With (11) and (7), we obtain that $\lim_{t \to 0} T(t, y; \psi) = \psi(y)$. With (9) and (10), $T(t, \cdot; \psi)$ satisfies the interface conditions (2) and (3) for any $t > 0$.

Thanks to the convenient transform (7), problem (8)-(11) is now formulated as a diffusive process in a discontinuous medium with generalized interface conditions as considered in [4]. Let us recall in the following box the results from [4] and useful in this context.
with an interface at 0 and an initial condition $\psi$:

$$ \phi(y)\partial_t T(t, y; \psi) = \mathcal{F}_y T(t, y; \psi), \quad (1) $$

$$ T(t, 0^{-}; \psi) = T(t, 0^{+}; \psi), \quad (2) $$

$$ k_1 \nabla T(t, 0^{-}; \psi) = k_2 \nabla T(t, 0^{+}; \psi), \quad (3) $$

$$ T(0, y; \psi) = \psi(y), \quad (4) $$

where $T(t, y; \psi)$ is the temperature at time $t$ and point $y$,

$$ \mathcal{F}_y f(y) := \nabla_y(k(y)\nabla_y f(y)) \quad \text{and} \quad \phi(y) := \rho(y)c(y). $$

Here, $\rho$, $c$ and $k$ are functions constant on each side of the interface located at 0. For such functions $f$, we write $f(y) = f_1$ if $y < 0$ and $f_2$ if $y > 0$. The constants are the thermal conductivity $k_i$ [W m$^{-1}$ K$^{-1}$], the specific heat $c_i$ [J kg$^{-1}$ K$^{-1}$], and the density $\rho_i$ [kg m$^{-3}$]. All these constants are assumed to be positive.

First, we consider the problem in free space, and later in a finite domain $[-L, L]$ with reflected boundary conditions at each side. The computations are the same in such a case.

Integrating (1) on the domain and using the Green formula,

$$ \begin{aligned} \partial_t \int_{-\infty}^{+\infty} \phi(y) T(t, y; \psi) \, dy &= \int_{-\infty}^{0} \nabla_y(k(y)\nabla_y T(t, y; \psi)) \, dy + \int_{0}^{+\infty} \nabla_y(k(y)\nabla_y T(t, y; \psi)) \, dy \\ &= k_1 \nabla T(t, 0^{-}; \psi) - k_2 \nabla T(t, 0^{+}; \psi) = 0. \quad (5) \end{aligned} $$

Hence, the quantity $\int_{-\infty}^{+\infty} \phi(y)T(t, y; \psi) \, dy$ is conserved over time. On the other hand we observe that,

$$ \begin{aligned} \partial_t \int_{-\infty}^{+\infty} T(t, y; \psi) \, dy &= \int_{-\infty}^{0} \frac{1}{\phi_1} \nabla_y(k_1 \nabla_y T(t, y; \psi)) \, dy + \int_{0}^{+\infty} \frac{1}{\phi_2} \nabla_y(k_2 \nabla_y T(t, y; \psi)) \, dy \\ &= \frac{k_1}{\phi_1} \nabla T(t, 0^{-}; \psi) - \frac{k_2}{\phi_2} \nabla T(t, 0^{+}; \psi). \quad (6) \end{aligned} $$

Therefore, the quantity $\int_{-\infty}^{+\infty} T(t, y; \psi) \, dy$ is not conserved over time.

## 2.2 Lagrangian approach

As practical implication of (5) we design a numerical scheme relying upon a Lagrangian approach where the thermal energy of the system is discretized by a set of particles, each one carrying a specific quantity. The design of such diffusive Lagrangian scheme is achieved by defining a fundamental solution $q(t, x, y)$ [m$^{-1}$] such that $q(t, x, y) \, dy$ is the probability that
In [4], we were considering solving the family of problems

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

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

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

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

with $D(y) = D_1$ if $y < 0$ and $D_2$ if $y > 0$. 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 (12)-(15) is called a fundamental solution.

A closed-form expression for $q$ is given by the method of images:

* For a negative source $(x < 0)$, the fundamental solution $q$ is

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

with $\beta := \frac{\sqrt{D_2}}{\sqrt{D_1}}$, $A \in (-1, 1)$ and $B > 0$, formulated as

$$ A = \frac{\Gamma - 1}{\Gamma + 1}, \quad \text{and} \quad B = \frac{1 - \gamma}{1 + \gamma} \frac{2\beta^2}{1 + \Gamma}, $$

where $\mu := \frac{\nu_2 - \nu_1}{\nu_2 + \nu_1}$, $\gamma := \frac{\kappa_2 - \kappa_1}{\kappa_2 + \kappa_1}$ and $\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}$.

* For a positive source $(x > 0)$, 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} \quad (17) $$

with $\tilde{\beta} := \frac{\sqrt{D_1}}{\sqrt{D_2}} = \frac{1}{\beta}$, $\tilde{A} \in (-1, 1)$ and $\tilde{B} > 0$, formulated as

$$ \tilde{A} = \frac{\tilde{\Gamma} - 1}{\tilde{\Gamma} + 1} = \frac{1 - \tilde{\Gamma}}{1 + \tilde{\Gamma}} = -A, \quad \text{and} \quad \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, $$

with $\mu := \frac{\nu_2 - \nu_1}{\nu_2 + \nu_1}$, $\gamma := \frac{\kappa_2 - \kappa_1}{\kappa_2 + \kappa_1}$ and $\tilde{\Gamma} := \tilde{\beta} \frac{1 + \gamma}{1 - \gamma} \frac{1 - \mu}{1 + \mu} = \frac{1}{\Gamma}$.

For a negative source (resp. for a positive source), the condition $A + B = 1$ (resp. $\tilde{A} + \tilde{B} = 1$) is equivalent to mass conservation. This condition is met if and only if

$$ \beta^2 = \frac{1}{\tilde{\beta}^2} = \frac{D_2}{D_1} = \frac{1 + \gamma}{1 - \gamma} = \frac{\kappa_2}{\kappa_1}. \quad (18) $$
a particle starting from $x$ ends in a small volume $dy$ around $y$ at time $t$. We relate $q$ and $T$ through

$$T(t, y; \psi) = \int \psi(x)q(t, x, y) \frac{\phi(x)}{\phi(y)} dx. \qquad (7)$$

Here, $q$ is the family of solutions to

$$\partial_t q(t, x, y) = \mathcal{L}_y q(t, x, y), \qquad (8)$$

$$\frac{1}{\phi_1} q(t, x, 0^{-}) = \frac{1}{\phi_2} q(t, x, 0^{+}), \qquad (9)$$

$$\alpha_1 \nabla q(t, x, 0^{-}) = \alpha_2 \nabla q(t, x, 0^{+}), \qquad (10)$$

$$\lim_{t \to 0} q(t, x, y) = \delta_x(y), \qquad (11)$$

where

$$\mathcal{L}_y f(y) = \nabla_y(\alpha(y)\nabla_y f(y)) \quad \text{with} \quad \alpha_i := \frac{k_i}{\rho_i c_i} [\text{W m}^2 \text{ J}^{-1}], \quad i = 1, 2,$$

for any suitable function $f$. The quantity $\alpha$ is the *thermal diffusivity*.

With (11) and (7), we obtain that $\lim_{t \to 0} T(t, y; \psi) = \psi(y)$. With (9) and (10), $T(t, \cdot; \psi)$ satisfies the interface conditions (2) and (3). We note that (8)-(11) is formulated as a diffusive process in a discontinuous medium with generalized interface conditions, as indicated in Section 2 of [2], by assuming

$$D_i = \alpha_i, \quad \nu_i = \frac{1}{\phi_i} = \frac{1}{\rho_i c_i} \quad \text{and} \quad \kappa_i = \alpha_i \text{ for } i = 1, 2.$$

Hence, a closed-form expression for $q$ is given by the method of images. The mass conservation property (12) holds.

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

We also note that $q(t, x, y)$ is positive for any $t > 0$ and any $x, y \in \mathbb{R}$ and satisfies both the Kolmogorov backward equation in (13) and the *Chapman-Kolmogorov equation*.

$$\partial_t q(t, x, y) = \nabla_x(\alpha(x)\nabla_x q(t, x, y)) \qquad (13)$$

with the interface conditions

$$q(t, 0^{-}, y) = q(t, 0^{+}, y) \quad \text{and} \quad k_1 \nabla_x q(t, 0^{-}, y) = k_2 \nabla_x q(t, 0^{+}, y) \quad \text{for any } y.$$
Now, by analogy between (8)-(11) and (12)-(15),

$$D_i = \alpha_i, \nu_i = \frac{1}{\phi_i} = \frac{1}{\rho_i c_i} \quad \text{and} \quad \kappa_i = \alpha_i \text{ for } i = 1, 2; \qquad (19)$$

The mass conservation property (18) holds.

We also note that $q(t, x, y)$ is positive for any $t > 0$ and any $x, y \in \mathbb{R}$ and satisfies both the Kolmogorov backward equation in (20) and the Chapman-Kolmogorov equation.

$$\partial_t q(t, x, y) = \nabla_x(\alpha(x)\nabla_x q(t, x, y)) \qquad (20)$$

with the interface conditions

$$q(t, 0-, y) = q(t, 0+, y) \quad \text{and} \quad k_1 \nabla_x q(t, 0-, y) = k_2 \nabla_x q(t, 0+, y) \quad \text{for any } y.$$

From the properties of the fundamental solution, $q(t, x, y)$ satisfies

$$q(t + s, x, y) = \int_{\mathbb{R}} q(s, x, z)q(t, z, y) \, dz \quad \text{for any } s, t > 0, x, y \in \mathbb{R}. \qquad (21)$$

The positivity and the mass conservation are sufficient to interpret $q(t, x, \cdot)$ as a probability density. The Chapman-Kolmogorov property ensures that moving the particle again on time steps $t$ and $s$ is the same as performing one step at time $s + t$ (the so-called *Markov property*).

Finally, since $\mathcal{L}_y\phi(y) = 0$, it follows from (20) that

$$\partial_t \int \phi(x)q(t, x, y) \, dx = \int \phi(x)\nabla_x(\alpha(x)\nabla_x q(t, x, y)) \, dx$$

and using the Green formula twice,

$$\partial_t \int \phi(x)q(t, x, y) \, dx = 0.$$

Hence, $\phi$ is an invariant measure for $q$. It corresponds to the *steady state*. With (7), $t \mapsto T(t, y; 1)$ is constant.

### 2.2.1 GUM, GSBM, and GHYMLA algorithms

In [4], three algorithms, namely GUM, GSBM, and GHYMLA algorithms, were proposed to simulate the density transition function $q$. Details on the theoretical derivation of these algorithms are provided in [4] and are reported in Supplementary Materials S1. Notice these algorithms are restricted to the one dimensional case, extending them to handle higher-dimensional cases is not straightforward.

Applied to heat transfer problem, using the analogy (19) and the definitions of $A, B$, (resp. $\tilde{A}, \tilde{B}$) provided in (17) (resp. (16)), the constants required by the computational algorithms (see Section S1 in the Supplementary Materials) are provided in Table 1, as a function of the physical properties of the involved materials.

These three algorithms allow to solve the same diffusive problem being based on different approaches that affect the behavior of the particle close to the discontinuity interface (interface layer).
From the properties of the fundamental solution, $q(t, x, y)$ satisfies

$$q(t + s, x, y) = \int_{\mathbb{R}} q(s, x, z)q(t, z, y) \, dz \text{ for any } s, t > 0, x, y \in \mathbb{R}. \quad (14)$$

The positivity and the mass conservation are sufficient to interpret $q(t, x, \cdot)$ as a probability density. The Chapman-Kolmogorov property ensures that moving the particle again on time steps $t$ and $s$ is the same as performing one step at time $s + t$ (the so-called *Markov property*).

Finally, since $\mathcal{L}_y\phi(y) = 0$, it follows from (13) that

$$\partial_t \int \phi(x)q(t, x, y) \, dx = \int \phi(x)\nabla_x(\alpha(x)\nabla_x q(t, x, y)) \, dx$$

and using the Green formula twice,

$$\partial_t \int \phi(x)q(t, x, y) \, dx = 0.$$

Hence, $\phi$ is an invariant measure for $q$. It corresponds to the *steady state*. With (7), $t \mapsto T(t, y; 1)$ is constant.

### 2.2.1 Setting up the Lagrangian method

To solve the one-dimensional thermal diffusion process, we design a Lagrangian method based on $N$ particles and applied to a finite domain $[-L, L]$ with reflection boundary conditions. As mentioned in [2], when we are restricted to the finite domain, the closed-form expression $q$ has to be replaced by a fundamental solution $q_L$ that still satisfies both the Kolmogorov forward equation and the Kolmogorov backward equation with the additional conditions

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

At the boundaries $\pm L$, $q(t, x, \cdot)$ approximates $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 constants required by the computational algorithms are provided in Table (1), as a function of the physical properties of the involved materials.

For each particle the position $X(0)$ is drawn from a prescribed given initial distribution $\nu$ over $[-L, L]$, in line with the considered initial condition. At any time $t$ the 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, \dots, n}$
<table><caption>Table 1: Constants in the algorithms for the Fourier's law with interface.</caption><thead><tr><th>Algorithm</th><th colspan="2">Parameters</th></tr></thead><tbody><tr><td rowspan="2">GUM</td><td>A = <div style="display: inline-block; vertical-align: middle; text-align: center;"><div style="border-bottom: 1px solid black; padding-bottom: 1px;">ρ<sub>1</sub>c<sub>1</sub>√α<sub>1</sub> - ρ<sub>2</sub>c<sub>2</sub>√α<sub>2</sub></div><div style="padding-top: 1px;">ρ<sub>1</sub>c<sub>1</sub>√α<sub>1</sub> + ρ<sub>2</sub>c<sub>2</sub>√α<sub>2</sub></div></div></td><td>B = <div style="display: inline-block; vertical-align: middle; text-align: center;"><div style="border-bottom: 1px solid black; padding-bottom: 1px;">2ρ<sub>2</sub>c<sub>2</sub>√α<sub>2</sub></div><div style="padding-top: 1px;">ρ<sub>1</sub>c<sub>1</sub>√α<sub>1</sub> + ρ<sub>2</sub>c<sub>2</sub>√α<sub>2</sub></div></div></td></tr><tr><td>˜A = <div style="display: inline-block; vertical-align: middle; text-align: center;"><div style="border-bottom: 1px solid black; padding-bottom: 1px;">ρ<sub>2</sub>c<sub>2</sub>√α<sub>2</sub> - ρ<sub>1</sub>c<sub>1</sub>√α<sub>1</sub></div><div style="padding-top: 1px;">ρ<sub>1</sub>c<sub>1</sub>√α<sub>1</sub> + ρ<sub>2</sub>c<sub>2</sub>√α<sub>2</sub></div></div></td><td>˜B = <div style="display: inline-block; vertical-align: middle; text-align: center;"><div style="border-bottom: 1px solid black; padding-bottom: 1px;">2ρ<sub>1</sub>c<sub>1</sub>√α<sub>1</sub></div><div style="padding-top: 1px;">ρ<sub>1</sub>c<sub>1</sub>√α<sub>1</sub> + ρ<sub>2</sub>c<sub>2</sub>√α<sub>2</sub></div></div></td></tr><tr><td>GSBM, GHYMLA</td><td colspan="2">θ<sub>g</sub> = <div style="display: inline-block; vertical-align: middle; text-align: center;"><div style="border-bottom: 1px solid black; padding-bottom: 1px;">ρ<sub>2</sub>c<sub>2</sub>√α<sub>2</sub> - ρ<sub>1</sub>c<sub>1</sub>√α<sub>1</sub></div><div style="padding-top: 1px;">ρ<sub>2</sub>c<sub>2</sub>√α<sub>2</sub> + ρ<sub>1</sub>c<sub>1</sub>√α<sub>1</sub></div></div></td></tr></tbody></table>

* The GUM algorithm is derived from the method of the images. Analytical expressions are employed to calculate the probability that a particle is reflected or crosses the discontinuity interface. The algorithm does not require the calculation of a hitting time but relies on an approximation of the probability distribution associated with a particle position.

* The GSBM describes the particle's motion around the discontinuity interface by a stochastic process grounded on the Skew Brownian motion [21]. The particle's position during each time increment is determined through two steps: (1) computation of the first hitting time, which is the time it takes for the particle to reach the interface, to determine whether the particle reaches the discontinuity during the time step and (2) the calculation of the particle's final position based on Step (1).

* The GHYMLA algorithm takes physics-based approach ("splitting up jump"), dividing the particle's jump into two components: (1) the displacement occurring within one layer until the discontinuity interface, and (2) the motion from discontinuity to the other layer.

Broadly speaking GHYMLA and GSBM both involve a two-step calculation, as opposed to GUM that solely relies on the sampling of an approximate particle location probability. The main difference between the GHYMLA and GSBM lies in their estimation of the first hitting time $τ$. In the GSBM $τ$ is determined according to the *ExactHittingTime* algorithm which is grounded on the Brownian Bridge concept which requires the execution of further steps compared to the GHYMLA algorithms, as displayed in Algorithm 2 in Supplementary Material. *ExactHittingTime* algorithm provides a more accurate quantification but at a higher computational cost. GHYMLA, on the other hand, uses a linear relationship between the particle's position at times $t$ and $t + Δt$ to estimate $τ$. This approach offers faster computations but lower accuracy [4].

### 2.2.2 Setting up the Lagrangian method

Now we have algorithms to deal with particle displacements around the discontinuity, we can set up a Lagrangian method.

To solve the one-dimensional thermal diffusion process, we design a Lagrangian method based on $N$ particles and applied to a finite domain $[-L, L]$ with reflection boundary conditions. As mentioned by [4], when we are restricted to the finite domain, the closed-form expression $q$ has to be replaced by a fundamental solution $q_L$ that still satisfies both the Kolmogorov forward equation and the Kolmogorov backward equation with the additional conditions

$$ \nabla_x q_L(t, \pm L, y) = \nabla_y q_L(t, x, \pm L) = 0. $$
<table><caption>Table 1: Constants in the algorithms for the Fourier's law with interface.</caption><thead><tr><th>Algorithm</th><th>Parameters</th></tr></thead><tbody><tr><td>GUM</td><td>A = <sup>ρ<sub>1</sub>c<sub>1</sub>√α<sub>1</sub>−ρ<sub>2</sub>c<sub>2</sub>√α<sub>2</sub></sup>⁄<sub>ρ<sub>1</sub>c<sub>1</sub>√α<sub>1</sub>+ρ<sub>2</sub>c<sub>2</sub>√α<sub>2</sub></sub> B = <sup>2ρ<sub>2</sub>c<sub>2</sub>√α<sub>2</sub></sup>⁄<sub>ρ<sub>1</sub>c<sub>1</sub>√α<sub>1</sub>+ρ<sub>2</sub>c<sub>2</sub>√α<sub>2</sub></sub><br>˜A = <sup>ρ<sub>2</sub>c<sub>2</sub>√α<sub>2</sub>−ρ<sub>1</sub>c<sub>1</sub>√α<sub>1</sub></sup>⁄<sub>ρ<sub>1</sub>c<sub>1</sub>√α<sub>1</sub>+ρ<sub>2</sub>c<sub>2</sub>√α<sub>2</sub></sub> ˜B = <sup>2ρ<sub>1</sub>c<sub>1</sub>√α<sub>1</sub></sup>⁄<sub>ρ<sub>1</sub>c<sub>1</sub>√α<sub>1</sub>+ρ<sub>2</sub>c<sub>2</sub>√α<sub>2</sub></sub></td></tr><tr><td>GSBM, GHYMLA</td><td>θ<sub>g</sub> = <sup>ρ<sub>2</sub>c<sub>2</sub>√α<sub>2</sub>−ρ<sub>1</sub>c<sub>1</sub>√α<sub>1</sub></sup>⁄<sub>ρ<sub>2</sub>c<sub>2</sub>√α<sub>2</sub>+ρ<sub>1</sub>c<sub>1</sub>√α<sub>1</sub></sub></td></tr></tbody></table>

performs a random walk, and the distribution of $X(t)$ is close to $q_L(t, X(0), \cdot)$. The particle's behavior is determined according to the region of the domain where the particle lies at each time step:

1. $X(t) \in I = [-\sqrt{6D_1\Delta t}, \sqrt{6D_2\Delta t}]$: the particle resides in the *interface layer* containing the interface point $X_I = 0$. The position $X(t+1)$ is sampled from $q(\Delta t, x, \cdot)$, or any of its approximations.

2. $X(t) \in B_{-L} = [-L, -L + \sqrt{6D_1\Delta t}]$ or $X(t) \in B_L = [L - \sqrt{6D_2\Delta t}, L]$: the particle belongs to the *boundary layer*. The next particle's position is sampled from a Gaussian step (GSBM) or a uniform step (GUM or GHYMLA). The particle is reflected when it reaches the boundary $x = \pm L$.

3. $X(t) \notin (I \cup B_{\pm L})$: $X(t+1)$ is determined by performing an usual step using a Gaussian (GSBM) and or a uniform (GUM or GHYMLA) step.

The size of the interfaces and boundary layers are function of $\Delta t$, the latter being selected so that $I \cap B_L = I \cap B_{-L} = \emptyset$.

Physical quantities of the domain can be recovered by properly counting the particles. We are given an initial profile $\psi$ which is an integrable, non-negative function on $[-L, L]$. For the $k$-th particle among the $N$ particles, we draw the initial position $X^{(k)}(0)$ of a particle according to the initial distribution $x \mapsto \phi(x)\psi(x)/E$ with $E := \int_{-L}^{L} \psi(x)\phi(x) dx$. The particle is then moved to $X^{(k)}(t)$ using one of scheme described above.

We associate with the $k$-th particle the “weight” $w_k := EN^{-1}\phi(X^{(k)}(t))^{-1}$.

Given $g: [-L, L] \to \mathbb{R}$, the mean operator is

$$ \mathcal{M}(g; t, N, \psi) := \sum_{k=1}^{N} w_k g(X^{(k)}(t)) = \frac{E}{N} \sum_{k=1}^{N} \frac{g(X^{(k)}(t))}{\phi(X^{(k)}(t))}. $$
At the boundaries $\pm L$, $q(t, x, \cdot)$ approximates $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.

For each particle the position $X(0)$ is drawn from a prescribed given initial distribution $\nu$ over $[-L, L]$, in line with the considered initial condition. At any time $t$ the 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 particle's behavior is determined according to the region of the domain where the particle lies at each time step:

1. $X(t) \in I := [-\sqrt{6D_1\Delta t}, \sqrt{6D_2\Delta t}]$: the particle resides in the *interface layer* containing the interface point $X_I = 0$. The position $X(t+1)$ is sampled from $q(\Delta t, x, \cdot)$, or any of its approximations.

2. $X(t) \in B_{-L} := [-L, -L + \sqrt{6D_1\Delta t}]$ or $X(t) \in B_L := [L - \sqrt{6D_2\Delta t}, L]$: the particle belongs to the *boundary layer*. The next particle's position is sampled from a Gaussian step (GSBM) or a uniform step (GUM or GHYMLA). The particle is reflected when it reaches the boundary $x = \pm L$.

3. $X(t) \notin (I \cup B_{\pm L})$: $X(t+1)$ is determined by performing a usual step using a Gaussian (GSBM) and or a uniform (GUM or GHYMLA) step.

The size of the interfaces and boundary layers are function of $\Delta t$, the latter being selected so that $I \cap B_L = I \cap B_{-L} = \emptyset$.

Physical quantities of the domain can be recovered by properly counting the particles. While in a conventional particle tracking approach each particle represents a mass, here it indicates a quantity of energy. By assigning a weight to each particle, the Lagrangian approach can be implemented to thermal systems. The weight allows to attribute to each particle an amount of energy based on the thermal properties of medium/layer and/or a possible thermal source within the system. We are given an initial profile $\psi$ which is an integrable, non-negative function on $[-L, L]$. For the $k$-th particle among the $N$ particles, we draw the initial position $X^{(k)}(0)$ of a particle according to the initial distribution $x \mapsto \phi(x)\psi(x)/E$ with $E := \int_{-L}^{L} \psi(x)\phi(x) dx$. The particle is then moved to $X^{(k)}(t)$ using one of scheme described above.

We associate with the $k$-th particle the “weight” $w_k := EN^{-1}\phi(X^{(k)}(t))^{-1}$.

Given $g: [-L, L] \to \mathbb{R}$, the *mean operator* is

$$ \mathcal{M}(g; t, N, \psi) := \sum_{k=1}^{N} w_k g(X^{(k)}(t)) = \frac{E}{N} \sum_{k=1}^{N} \frac{g(X^{(k)}(t))}{\phi(X^{(k)}(t))}. $$

According to the strong law of large numbers (and assuming that we use an exact scheme),

$$ \begin{align} \lim_{N \to \infty} \mathcal{M}(g; t, N, \psi) &= E \int_{-L}^{L} \int_{-L}^{L} \frac{\phi(x)\psi(x)}{E} q_L(t, x, y) \frac{g(y)}{\phi(y)} dy dx \nonumber \\ &= \int_{-L}^{L} \int_{-L}^{L} \psi(x) q_L(t, x, y) \frac{\phi(x)}{\phi(y)} g(y) dy dx \tag{22} \\ &= \int_{-L}^{L} Tg(y) dy. \nonumber \end{align} $$
According to the strong law of large numbers (and assuming that we use an exact scheme),

$$ 
\begin{aligned}
\lim_{N \to \infty} \mathcal{M}(g; t, N, \psi) &= E \int_{-L}^{L} \int_{-L}^{L} \frac{\phi(x)\psi(x)}{E} q_L(t, x, y) \frac{g(y)}{\phi(y)} \, dy \, dx \\
&= \int_{-L}^{L} \int_{-L}^{L} \psi(x) q_L(t, x, y) \frac{\phi(x)}{\phi(y)} g(y) \, dy \, dx \quad (15) \\
&= \int_{-L}^{L} T(t, y; \psi) g(y) \, dy.
\end{aligned}
$$

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$.

We take $g_i(y) := \frac{1}{|A_i|} \mathbb{1}_{A_i}(y)$. Then

$$ \mathcal{M}(g_i; t, N, \psi) \approx \frac{1}{|A_i|} \frac{E}{\phi(a_i)} \frac{n(t, A_i)}{N} \quad (16) $$

with

$$ n(t, A) := \sum_{i=1}^{N} \mathbb{1}_{A}(X^{(k)}(t)) \quad (17) $$

# 3 Benchmark tests

Heat transfer in a one-dimensional bimaterial medium of finite size $[-L, L]$ with reflecting boundary conditions (RBC) as displayed in Figure 1 is simulated by the GUM Algorithm (Algorithm 1 in Supplementary Material), GSBM Algorithm (Algorithm 2 in Supplementary Material) and GHYMLA Algorithm (Algorithm 3 in Supplementary Material). 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, such as density ($\rho$), thermal conductivity ($k$), specific heat ($c$), and thermal diffusivity ($\alpha$), are constant within $R_1$ and $R_2$. The algorithms are tested by considering two case studies:

## Case 1: Thermal equilibrium

*Domain*: The two layers are initialized in thermal equilibrium, i.e., (20)-(21) holds.

*Aim*: Check the conservation of the equilibrium temperature of each layer and the thermal energy.

## Case 2: Thermal gradient

*Domain*: An initial temperature difference is applied only on the region $R_2$.
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$.

We take $g_i(y) := \frac{1}{|A_i|}1_{A_i}(y)$. Then

$$ \mathcal{M}(g_i; t, N, \psi) \approx \frac{1}{|A_i|} \frac{E}{\phi(a_i)} \frac{n(t, A_i)}{N} \qquad (23) $$

with

$$ n(t, A) := \sum_{i=1}^{N} 1_{A}(X^{(k)}(t)). \qquad (24) $$

# 3 Benchmark tests

Heat transfer in a one-dimensional bimaterial medium of finite size $[-L, L]$ with reflecting boundary conditions (RBC) as displayed in Figure 1 is simulated by the GUM Algorithm (Algorithm 1 in Supplementary Material), GSBM Algorithm (Algorithm 2 in Supplementary Material) and GHYMLA Algorithm (Algorithm 3 in Supplementary Material). 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, such as density ($\rho$), thermal conductivity ($k$), specific heat ($c$), and thermal diffusivity ($\alpha$), are constant within $R_1$ and $R_2$. The algorithms are tested by considering two case studies:

## Case 1: Thermal equilibrium

*Domain*: The two layers are initialized in thermal equilibrium, i.e., (28)-(29) holds.

*Aim*: Check the conservation of the equilibrium temperature of each layer and the thermal energy.

## Case 2: Thermal gradient

*Domain*: An initial temperature difference is applied only on the region $R_2$.

*Aim*: Check that the temperature distribution along the domain is consistent with the analytical solution (31) over time.

All benchmark tests have the same domain size and physical properties as those described in Table 2, but their initial conditions differ. The model parameters chosen for this study reflect typical values observed for conductive heat transfer at the water-soil interface, enabling testing the algorithms' performance in a realistic scenario. Each test assumes the domain boundaries at $x = \pm L$ are reflective barriers (i.e., no-flux boundary conditions). The algorithms for both test cases are publicly available in the repository [3].

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 are indicated in Table 3.

## 3.1 Case 1: Thermal equilibrium

The test is performed to verify that the GUM, GSBM and GHYMLA algorithms maintain the temperature within each layer constant over time under steady-state equilibrium condition. Here we provide the theoretical benchmark results to assess the Lagrangian formulations outlined in Section 2.2 followed by the numerical results.
Aim: Check that the temperature distribution along the domain is consistent with the analytical solution (23) over time.

All benchmark tests have the same domain size and physical properties as those described in Table 2, but their initial conditions differ. Each test assumes the domain boundaries at $x = \pm L$ are reflective barriers (i.e., no-flux boundary conditions).

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 are indicated in Table 3.

<table><caption>Table 2: Size and physical properties of the domain used in the tests.</caption><thead><tr><th></th><th><i>R</i><sub>1</sub></th><th><i>R</i><sub>2</sub></th><th>parameter</th><th>SI units</th></tr></thead><tbody><tr><td><i>k</i></td><td>0.603</td><td>0.335</td><td>thermal conductivity</td><td>W m<sup>-1</sup> K<sup>-1</sup></td></tr><tr><td><i>ρ</i></td><td>1000</td><td>1600</td><td>density</td><td>kg m<sup>-3</sup></td></tr><tr><td><i>c</i></td><td>4184</td><td>753</td><td>specific heat</td><td>J kg<sup>-1</sup> K<sup>-1</sup></td></tr><tr><td><i>L</i></td><td>0.005</td><td>0.005</td><td>length of the domain region</td><td>m</td></tr><tr><td><i>X</i><sub>I</sub></td><td>0</td><td></td><td>position of the interface</td><td>m</td></tr></tbody></table>

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

## 3.1 Case 1: Thermal equilibrium

The test is performed to verify that the GUM, GSBM and GHYMLA algorithms maintain the temperature within each layer constant over time under steady-state equilibrium condition. Here we provide the theoretical benchmark results to assess the Lagrangian formulations outlined in 2.2 followed by the numerical results.
<table><caption>Table 2: Size and physical properties of the domain used in the tests.</caption><thead><tr><th></th><th>R<sub>1</sub></th><th>R<sub>2</sub></th><th>parameter</th><th>SI units</th></tr></thead><tbody><tr><td>k</td><td>0.603</td><td>0.335</td><td>thermal conductivity</td><td>W m<sup>-1</sup> K<sup>-1</sup></td></tr><tr><td>ρ</td><td>1000</td><td>1600</td><td>density</td><td>kg m<sup>-3</sup></td></tr><tr><td>c</td><td>4184</td><td>753</td><td>specific heat</td><td>J kg<sup>-1</sup> K<sup>-1</sup></td></tr><tr><td>L</td><td>0.005</td><td>0.005</td><td>length of the domain region</td><td>m</td></tr><tr><td>X<sub>I</sub></td><td>0</td><td></td><td>position of the interface</td><td>m</td></tr></tbody></table>

<table><caption>Table 3: 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>

### 3.1.1 Theoretical benchmark solution

We deduce from the Kolmogorov backward equation presented by [4] that

$$ \partial_t T(t, y; \psi) = \partial_t \int_{-L}^{L} \psi(x) q_L(t, x, y) \frac{\phi(x)}{\phi(y)} dx = \int_{-L}^{L} \psi(x) \frac{\phi(x)}{\phi(y)} \nabla_x (\alpha(x) \nabla_x q_L(t, x, y)) dx. $$

Assume that $\psi(x) = \psi_1$ if $x < 0$ and $\psi(x) = \psi_2$ if $x > 0$. From the Green formula,

$$ \begin{aligned} & \int_{-L}^{L} \psi(x) \frac{\phi(x)}{\phi(y)} \nabla_x (\alpha(x) \nabla_x q_L(t, x, y)) dx \\ &= \int_{-L}^{0} \psi(x) \frac{\phi(x)}{\phi(y)} \nabla_x (\alpha(x) \nabla_x q_L(t, x, y)) dx + \int_{0}^{L} \psi(x) \frac{\phi(x)}{\phi(y)} \nabla_x (\alpha(x) \nabla_x q_L(t, x, y)) dx \\ &= \frac{1}{\phi(y)} (\psi_1 \phi_1 \alpha_1 \nabla_x q_L(t, 0-, y) - \psi_2 \phi_2 \alpha_2 \nabla_x q_L(t, 0+, y)), \quad (25) \end{aligned} $$

since $\nabla_x q_L(t, \pm L, y) = 0$ as we use the reflected boundary conditions. Owing to the interface conditions of the Kolmogorov backward equation,

$$ \begin{aligned} \partial_t T &= \frac{1}{\phi(y)} (\psi_1 \phi_1 \alpha_1 \nabla_x q_L(t, 0-, y) - \psi_2 \phi_2 \alpha_2 \nabla_x q_L(t, 0+, y)) \\ &= \frac{1}{\phi(y)} (\psi_1 - \psi_2) k_2 \nabla_x q_L(t, 0+, y). \quad (26) \end{aligned} $$

We then obtain that

$$ T = \psi(y) + (\psi_1 - \psi_2) \frac{k_2}{\phi(y)} \int_{0}^{t} \nabla_x q_L(s, 0+, y) \mathrm{d}s. \quad (27) $$

When $\psi$ equal to a constant $T_0$, $T(t, y; T_0) = T_0$ for any $y$, being also an immediate consequence from (27). This implies that

$$ \int_{-L}^{0} T(t, y; T_0) \mathrm{d}y = \int_{-L}^{0} T(0, y; T_0) \mathrm{d}y = LT_0, \quad (28) $$

$$ \int_{L}^{1} T(t, y; T_0) \mathrm{d}y = \int_{0}^{L} T(0, y; T_0) \mathrm{d}y = LT_0. \quad (29) $$
### 3.1.1 Theoretical benchmark solution

We deduce from the Kolmogorov backward equation presented in [2] that

$$ \partial_t T(t, y; \psi) = \partial_t \int_{-L}^{L} \psi(x) q_L(t, x, y) \frac{\phi(x)}{\phi(y)} dx = \int_{-L}^{L} \psi(x) \frac{\phi(x)}{\phi(y)} \nabla_x (\alpha(x) \nabla_x q_L(t, x, y)) dx. $$

Assume that $\psi(x) = \psi_1$ if $x < 0$ and $\psi(x) = \psi_2$ if $x > 0$. From the Green formula,

$$ \begin{aligned} & \int_{-L}^{L} \psi(x) \frac{\phi(x)}{\phi(y)} \nabla_x (\alpha(x) \nabla_x q_L(t, x, y)) dx \\ &= \int_{-L}^{0} \psi(x) \frac{\phi(x)}{\phi(y)} \nabla_x (\alpha(x) \nabla_x q_L(t, x, y)) dx + \int_{0}^{L} \psi(x) \frac{\phi(x)}{\phi(y)} \nabla_x (\alpha(x) \nabla_x q_L(t, x, y)) dx \\ &= \frac{1}{\phi(y)} (\psi_1 \phi_1 \alpha_1 \nabla_x q_L(t, 0-, y) - \psi_2 \phi_2 \alpha_2 \nabla_x q_L(t, 0+, y)), \quad (18) \end{aligned} $$

since $\nabla_x q_L(t, \pm L, y) = 0$ as we use the reflected boundary conditions. Owing to the interface conditions of the Kolmogorov backward equation,

$$ \begin{aligned} \partial_t T(t, y; \psi) &= \frac{1}{\phi(y)} (\psi_1 \phi_1 \alpha_1 \nabla_x q_L(t, 0-, y) - \psi_2 \phi_2 \alpha_2 \nabla_x q_L(t, 0+, y)) \\ &= \frac{1}{\phi(y)} (\psi_1 - \psi_2) k_2 \nabla_x q_L(t, 0+, y). \end{aligned} $$

We then obtain that

$$ T(t, y; \psi) = \psi(y) + (\psi_1 - \psi_2) \frac{k_2}{\phi(y)} \int_{0}^{t} \nabla_x q_L(s, 0+, y) \, ds. \quad (19) $$

When $\psi$ equal to a constant $T_0$, $T(t, y; T_0) = T_0$ for any $y$, being also an immediate consequence from (19). This implies that

$$ \int_{-L}^{0} T(t, y; T_0) \, dy = \int_{-L}^{0} T(0, y; T_0) \, dy = LT_0, \quad (20) $$

$$ \int_{L}^{1} T(t, y; T_0) \, dy = \int_{0}^{L} T(0, y; T_0) \, dy = LT_0. \quad (21) $$

Using the notation (17) with $R_1 = [-L, 0]$ and $R_2 = [0, L]$, we estimate these quantities by

$$ \begin{aligned} \int_{0}^{L} T(t, y; T_0) \, dy &\approx \mathcal{M}(1_{[0,L]}; t, N, T_0) = T_0 L \frac{\phi_1 + \phi_2 n(t, R_2)}{\phi_2 N}, \\ \int_{-L}^{0} T(t, y; T_0) \, dy &\approx \mathcal{M}(1_{[-L,0]}; t, N, T_0) = T_0 L \frac{\phi_1 + \phi_2 n(t, R_1)}{\phi_1 N} \end{aligned} $$
Using the notation (24) with $R_1 = [-L, 0]$ and $R_2 = [0, L]$, we estimate these quantities by

$$ \int_{0}^{L} T(t, y; T_0) \, dy \approx \mathcal{M}(1_{[0,L]}; t, N, T_0) = T_0 L \frac{\phi_1 + \phi_2 n(t, R_2)}{\phi_2 N}, \\ \int_{-L}^{0} T(t, y; T_0) \, dy \approx \mathcal{M}(1_{[-L,0]}; t, N, T_0) = T_0 L \frac{\phi_1 + \phi_2 n(t, R_1)}{\phi_1 N} $$

when $X(0)$ is distributed using the measure $\phi(x)/\int_{-L}^{L} \phi(x) \, dx$. Since $\phi$ is piecewise constant, the distributions of $X(0)$ given $X(0) < 0$ is uniform over $[-L, 0]$ and the distribution of $X(0)$ given $X(0) > 0$ is uniform over $[0, L]$. Besides,

$$ \mathbb{P}[X(0) \in R_i] = \frac{\phi_i}{\phi_1 + \phi_2} \approx \frac{n(0, R_i)}{N} \quad \text{for } i = 1, 2. $$

Therefore,

$$ \frac{n(0, R_1)}{n(0, R_2)} \approx \frac{\phi_1}{\phi_2} = \frac{\rho_1 c_1}{\rho_2 c_2}. $$

To check (28)-(29), it is sufficient to check that

$$ \frac{n(t, R_1)}{n(0, R_1)} \approx 1 \quad \text{and} \quad \frac{n(t, R_2)}{n(0, R_2)} \approx 1 \quad \text{for any } t \ge 0. \qquad (30) $$

### 3.1.2 Numerical results

The population of $N$ particles is split in two sub-populations with respective sizes $N_1$ and $N_2$ such that $N_1 + N_2 = N$ and $N_1/N_2 = \rho_1 c_1/\rho_2 c_2$. The particles from the first sub-population are uniformly distributed on $[-L, 0]$, the ones of the second sub-population are uniformly distributed on $[0, L]$, as displayed in Figure 1. We then check the conservation of the thermal equilibrium by using (30).

<figure><img src="image_1.png" /><figcaption>Figure 1: Case 1 setup: blue and red-filled regions indicate $R_1, R_2$, the continuous black lines indicate initial uniform particle placement in the two layers, different line thickness is used to indicate different particle density.</figcaption></figure>

In our numerical test we set $N_1 = 2 \times 10^6$ and $N_2 = 5.759 \times 10^5$. For all algorithms, simulations are run until $t = 10$ s with the same time step equal to $10^{-2}$ s to facilitate comparison.

The standard deviation $\sigma$ of the ratio $n(t, R_i)/N_i$ is used as indicator to evaluate the performance of the algorithm, lowest values of $\sigma$ indicating better performance. Results in Table 4 indicate that the three algorithms provide a standard deviation of the same order of magnitude with the highest value observed in the GHYMLA. We observe no significant trends in the results reported in Figure 2 thus we conclude that the GUM, GSBM and GHYMLA algorithms guarantee thermal energy conservation and maintain the temperature constant under steady-state equilibrium conditions. Additional analysis in Supplementary Material S2 reveals no evidence of bias in our numerical results.
when $X(0)$ is distributed using the measure $\phi(x)/\int_{-L}^{L} \phi(x) dx$. Since $\phi$ is piecewise constant, the distributions of $X(0)$ given $X(0) < 0$ is uniform over $[-L, 0]$ and the distribution of $X(0)$ given $X(0) > 0$ is uniform over $[0, L]$. Besides,

$$ \mathbb{P}[X(0) \in R_i] = \frac{\phi_i}{\phi_1 + \phi_2} \approx \frac{n(0, R_i)}{N} \quad \text{for } i = 1, 2. $$

Therefore,

$$ \frac{n(0, R_1)}{n(0, R_2)} \approx \frac{\phi_1}{\phi_2} = \frac{\rho_1 c_1}{\rho_2 c_2}. $$

To check (20)-(21), it is sufficient to check that

$$ \frac{n(t, R_1)}{n(0, R_1)} \approx 1 \quad \text{and} \quad \frac{n(t, R_2)}{n(0, R_2)} \approx 1 \quad \text{for any } t \ge 0. \qquad (22) $$

### 3.1.2 Numerical results

The population of $N$ particles is split in two sub-populations with respective sizes $N_1$ and $N_2$ such that $N_1 + N_2 = N$ and $N_1/N_2 = \rho_1 c_1 / \rho_2 c_2$. The particles from the first sub-population are uniformly distributed on $[-L, 0]$, the ones of the second sub-population are uniformly distributed on $[0, L]$, as displayed in Figure 1. We then check the conservation of the thermal equilibrium by using (22).

<figure><img src="image_3.png" /><figcaption>Figure 1: Case 1 setup: blue and red-filled regions indicate $R_1, R_2$, the continuous black lines indicate initial uniform particle placement in the two layers, different line thickness is used to indicate different particle density.</figcaption></figure>

In our numerical test we set $N_1 = 2 \times 10^6$ and $N_2 = 5.759 \times 10^5$. For all algorithms, simulations are run until $t = 1$ s with the same time step equal to $10^{-2}$ s to facilitate comparison.

The standard deviation $\sigma$ of the ratio $n(t, R_i)/N_i$ is used as indicator to evaluate the performance of the algorithm, lowest values of $\sigma$ indicating better performance. Results in Table 4 indicate that the three algorithms provide a standard deviation of the same order of magnitude with the highest value observed in the GHYMLA. We observe no significant trends in the results reported in Figure 2 thus we conclude that the GUM, GSBM and GHYMLA algo-
<figure><img src="image_2.png" /></figure>

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

<figure><img src="image_4.png" /><figcaption>Figure 2: Case 1: Temporal variation of the ratio $n(t, R_i)/N_i$ between the number of particles at time $t$ and $t = 0$ in $R_1$ (blue line) and $R_2$ (red line) for (a) GUM, (b) GSBM and (c) GHYMLA algorithm.</figcaption></figure>
rithms guarantee thermal energy conservation and maintain the temperature constant under steady-state equilibrium conditions.

<table><caption>Table 4: Case 1: Standard deviation $\sigma$ of $n(t, R_i)/N_i$ for the GUM, GSBM and GHYMLA algorithm.</caption><thead><tr><th></th><th>GUM</th><th>GSBM</th><th>GHYMLA</th></tr></thead><tbody><tr><td>σ</td><td>2.342 × 10<sup>-4</sup></td><td>2.316 × 10<sup>-4</sup></td><td>2.681 × 10<sup>-4</sup></td></tr></tbody></table>

## 3.2 Case 2: Thermal gradient

The test case simulates heat transport and the local temperature evolution in a bimaterial medium subjected to an initial temperature difference (see Figure 3) using the GUM, GSBM and GHYMLA algorithms. The numerical results are presented here after a theoretical formulation of the physical problem.

### 3.2.1 Theoretical benchmark solution

We follow the approach specified below. Set $\psi(x) = T_01_{[0,L]}(x)$, that is the Heaviside function, so that $\psi_1 = 0$, $\psi_2 = T_0$ for $T_0$ a temperature in K. Replacing $q_L$ by $q$ in (19),

$$T(t, y; T_0 1_{[0,L]}) = T_0 1_{[0,L]}(y) - \frac{k_2}{\phi(y)} T_0 \int_0^t \nabla_x q_L(s, 0+, y) \, ds.$$

Replacing $q_L(s, 0+, y)$ by $q(s, 0+, y)$ and using the closed-form of the density transition function $q$ in [2] we obtain

$$T(t, y; T_0 1_{[0,L]}) \approx T_0 \frac{\sqrt{k_2 \rho_2 c_2}}{\sqrt{k_2 \rho_2 c_2 + \sqrt{k_1 \rho_1 c_1}}} \times \begin{cases} \operatorname{erfc}\left(\frac{|y|}{2\sqrt{\alpha_1 t}}\right) & \text{for } y < 0, \\ 1 + \frac{\sqrt{k_1 \rho_1 c_1}}{\sqrt{k_2 \rho_2 c_2}} \operatorname{erf}\left(\frac{y}{2\sqrt{\alpha_2 t}}\right) & \text{for } y > 0. \end{cases} \quad (23)$$

Such an expression is already found in [3].

On the other hand, using the notations of Section 2.2.1,

$$T(t, a_i; T_0 1_{[0,T]}) \approx \int_{-L}^{L} T(t, y; T_0 1_{[0,L]}) g_i(y) \, dy \\ \approx \mathcal{M}(g_i; t, N, T_0 1_{[0,L]}) = T_0 \frac{L}{|A_i|} \frac{\phi_2}{\phi(a_i)} \frac{n(t, A_i)}{N}. \quad (24)$$
<table><caption>Table 4: Case 1: Standard deviation σ of $n(t, R_i)/N_i$ for the GUM, GSBM and GHYMLA algorithm.</caption><thead><tr><th></th><th>GUM</th><th>GSBM</th><th>GHYMLA</th></tr></thead><tbody><tr><th>σ</th><td>5.45 × 10<sup>-4</sup></td><td>5.36 × 10<sup>-4</sup></td><td>6.12 × 10<sup>-4</sup></td></tr></tbody></table>

## 3.2 Case 2: Thermal gradient

The test case simulates heat transport and the local temperature evolution in a bimaterial medium subjected to an initial temperature difference (see Figure 3) using the GUM, GSBM and GHYMLA algorithms. The numerical results are presented here after a theoretical formulation of the physical problem.

### 3.2.1 Theoretical benchmark solution

We follow the approach specified below. Set $\psi(x) = T_0 1_{[0,L]}(x)$, that is the Heaviside function, so that $\psi_1 = 0$, $\psi_2 = T_0$ for $T_0$ a temperature in K. Replacing $q_L$ by $q$ in (27),

$$T(t, y; T_0 1_{[0,L]}) = T_0 1_{[0,L]}(y) - \frac{k_2}{\phi(y)} T_0 \int_0^t \nabla_x q_L(s, 0+, y) \, ds.$$

Replacing $q_L(s, 0+, y)$ by $q(s, 0+, y)$ and using the closed-form of the density transition function $q$ presented by [4] we obtain

$$T(t, y; T_0 1_{[0,L]}) \approx T_0 \frac{\sqrt{k_2 \rho_2 c_2}}{\sqrt{k_2 \rho_2 c_2} + \sqrt{k_1 \rho_1 c_1}} \times \begin{cases} \operatorname{erfc}\left(\frac{|y|}{2\sqrt{\alpha_1 t}}\right) & \text{for } y < 0, \\ 1 + \frac{\sqrt{k_1 \rho_1 c_1}}{\sqrt{k_2 \rho_2 c_2}} \operatorname{erf}\left(\frac{y}{2\sqrt{\alpha_2 t}}\right) & \text{for } y > 0. \end{cases} \quad (31)$$

Such an expression was already proposed by [5].

On the other hand, using the notations of Section 2.2.2,

$$T(t, a_i; T_0 1_{[0,T]}) \approx \int_{-L}^{L} T(t, y; T_0 1_{[0,L]}) g_i(y) \, dy \\ \approx \mathcal{M}(g_i; t, N, T_0 1_{[0,L]}) = T_0 \frac{L}{|A_i|} \frac{\phi_2}{\phi(a_i)} \frac{n(t, A_i)}{N}. \quad (32)$$

### 3.2.2 Numerical results

The temperature is initially assumed step-wise constant and equal to 0 K in $R_1$ and $T_0 = 283.15$ K in $R_2$. The particles are initially uniformly distributed over $R_2$, as indicated in Figure 3. The reference temperature $T_{\text{ref}}(t, a_i)$ is computed according to (31) with $y = a_i$, while the numerical temperature is defined as

$$T_{\text{num}}(t, a_i) := T_0 \frac{L}{|A_i|} \frac{\phi_2}{\phi(a_i)} \frac{n(t, A_i)}{N}. \quad (33)$$
<figure><img src="image_4.png" /><figcaption>Figure 2: Case 1: Temporal variation of the ratio $n(t, R_i)/N_i$ between the number of particles at time $t$ and $t = 0$ in $R_1$ (blue line) and $R_2$ (red line) for (a) GUM, (b) GSBM and (c) GHYMLA algorithm.</figcaption></figure>
The quality of the solution is then evaluated through the root mean squared error applied to the temperature spatial distribution as different time levels

$$ \text{RMSE} = \sqrt{\frac{\sum_{i=1}^{N_a} (T_{\text{ref}}(t, a_i) - T_{\text{num}}(t, a_i))^2}{N_a}} \quad \text{in K.} \qquad (34) $$

<figure><img src="image_5.png" /><figcaption>Figure 3: Case 2 setup: blue and red regions indicate $R_1$, $R_2$, the continuous black line indicates initial uniform particle placement in layer $R_2$.</figcaption></figure>

The accuracy of temperature estimations across the domain is strongly affected by both the time increment magnitude and the width of the interface layer, the latter being a function of the chosen time step and the diffusive properties of the domain. Additionally, the proposed computational approach, relying upon the particle tracking method, demonstrates an inherent dependence on the number of particles employed in the simulation. The influence of the time step on the quantification of temperature over the domain is evaluated by the RMSE between $T_{\text{ref}}$ and $T_{\text{num}}$ for different values of the time step $\Delta t$. As shown in Figure 4, the GSBM accuracy is essentially independent of $\Delta t$ over the span of three orders of magnitudes. The trend between $\Delta t$ and RMSE observed in Figure 4 for both GHYMLA and GUM algorithms indicates a sensitivity of the algorithms to the chosen time step. The GUM algorithm attains a stable value of the error for $\Delta t \le 0.1$ s, while a sharp nonlinear increase of the error is observed for $\Delta t > 0.1$ s. This result suggests the existence of a threshold value for which we observe a correlation between RMSE and $\Delta t$ for GUM, in line with [4]. However, the threshold value of $\Delta t$ for which the error starts increasing is here different, suggesting that the threshold may be case-dependent. A direct linear convergence is observed for GHYMLA, i.e., the RMSE decreases linearly with $\Delta t$ like $\text{RMSE} \propto \Delta t^{0.4}$. For all tested values of $\Delta t$ the GHYMLA algorithm is the least accurate.

<figure><img src="image_6.png" /><figcaption>Figure 4: Case 2: RMSE between the analytical solution $T_{\text{ref}}$ and the numerical $T_{\text{num}}$ obtained from the GUM, GSBM and GHYMLA for different $\Delta t$ at time $t = 1$ s.</figcaption></figure>
### 3.2.2 Numerical results

The temperature is initially assumed step-wise constant and equal to 0 K in $R_1$ and $T_0 = 283.15$ K in $R_2$. The particles are initially uniformly distributed over $R_2$, as indicated in Figure 3. The reference temperature $T_{ref}(t, a_i)$ is computed according to (23) with $y = a_i$, while the numerical temperature is defined as

$$T_{\text{num}}(t, a_i) := T_0 \frac{L}{|A_i|} \frac{\phi_2}{\phi(a_i)} \frac{n(t, A_i)}{N}. \qquad (25)$$

The quality of the solution is then evaluated through the root mean squared error applied to the temperature spatial distribution as different time levels

$$\text{RMSE} = \sqrt{\frac{\sum_{i=1}^{N_a} (T_{\text{ref}}(t, a_i) - T_{\text{num}}(t, a_i))^2}{N_a}} \quad \text{in K.} \qquad (26)$$

<figure><img src="image_5.png" /><figcaption>Figure 3: Case 2 setup: blue and red regions indicate $R_1$, $R_2$, the continuous black line indicates initial uniform particle placement in layer $R_2$.</figcaption></figure>

The magnitude of the time step strongly affects accuracy. The influence of the time step on the quantification of temperature over the domain is evaluated by the RMSE between $T_{ref}$ and $T_{num}$ for different values of the time step $\Delta t$. As shown in Figure 4, the GSBM accuracy is essentially independent of $\Delta t$ over the span of three orders of magnitudes. The GUM algorithm attains a stable value of the error for $\Delta t \le 0.1$ s, while a sharp nonlinear increase of the error is observed for $\Delta t > 0.1$ s. This result suggests the existence of a threshold value for which we observe a correlation between RMSE and $\Delta t$ for GUM, in line with [2]. However, the threshold value of $\Delta t$ for which the error starts increasing is here different, suggesting that the threshold may be case-dependent.

A direct linear convergence is observed for GHYMLA, i.e., the RMSE decreases linearly with $\Delta t$. For all tested values of $\Delta t$ the GHYMLA algorithm is the least accurate.
The comparison between the analytical and the numerical temperature is displayed in Figure 5. The simulations are performed using the computational parameters in Table 3. The RMSE values at time $t = 2$ s obtained from the GUM, GSBM and GHYMLA are reported in Table 5 along with the average computational time $t_{cm}$. GSBM is the most accurate method, but has additional computational costs. GUM provides an approximated version of the transition density that speeds up the computation while maintaining good accuracy. GHYMLA has computational cost comparable to GUM, but is less accurate. Our findings in line with [4] support the conclusion that the GUM offers an effective trade-off between accuracy and computational efficiency.

<table><caption>Table 5: RMSE between the analytical and the numerical temperature at time $t = 2$ s and the average computational time $t_{cm}$ [s]</caption><thead><tr><th></th><th>RMSE</th><th>t<sub>cm</sub></th></tr></thead><tbody><tr><td>GUM</td><td>1.06</td><td>0.110</td></tr><tr><td>GSBM</td><td>0.849</td><td>12.087</td></tr><tr><td>GHYMLA</td><td>2.124</td><td>0.100</td></tr></tbody></table>
<figure><img src="image_6.png" /><figcaption>Figure 4: Case 2: RMSE between the analytical solution $T_{\text{ref}}$ and the numerical $T_{\text{num}}$ obtained from the GUM (red circles), GSBM (blue squares) and GHYMLA (black diamonds) for different $\Delta t$ at time $t = 1$ s.</figcaption></figure>

The comparison between the analytical and the numerical temperature is displayed in Figure 5. The simulations are performed using the computational parameters in Table 3.
<figure><img src="image_7.png" /></figure>

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

<figure><img src="image_9.png" /><figcaption>Figure 5: Case 2: Comparison between the analytical $T_{\text{ref}}$ (black dashed line) and the numerical $T_{\text{num}}$ (red line) temperature obtained from (a) GUM, (b) GSBM and (c) GHYMLA algorithm from the initial time to $t = 5$ s with a time interval of 1 s.</figcaption></figure>

A good accuracy in simulating the analytical solution is observed for all the assessed algorithms. The RMSE between $T_{\text{ref}}$ and $T_{\text{num}}$ is computed for an increasing number of particles $N_p$ as displayed in Figure 6. The GUM, GSBM, GHYMLA algorithms converge to the line RMSE $\propto N_p^{-0.5}$, i.e. the expected Monte Carlo convergence rate.
<figure><img src="image_7.png" /><figcaption>Figure 5: Case 2: Comparison between the analytical $T_{\text{ref}}$ (black dashed line) and the numerical $T_{\text{num}}$ (red line) temperature obtained from (a) GUM, (b) GSBM and (c) GHYMLA algorithm from the initial time to $t = 5$ s with a time interval of 1 s.</figcaption></figure>
<figure><img src="image_10.png" /><figcaption>Figure 6: Case 2: RMSE between the analytical $T_{\text{ref}}$ and the numerical $T_{\text{num}}$ temperature obtained from GUM, GSBM and GHYMLA for an increasing number of particles at time $t = 1$ s. The dashed line corresponds to RMSE $\propto N_p^{-0.5}$.</figcaption></figure>

Finally, we test the performance of the algorithms in estimating the temperature for different ratios between the thermal diffusivities of the two materials. The same thermal scenario described in Case 2 is simulated varying the thermal diffusivity gradient ($\alpha_1/\alpha_2$) over a two-orders-of-magnitude difference. The RMSE between the analytical and numerical temperature is used as indicator to estimate the accuracy of the algorithms. As displayed in Figure 7, the RMSE for each algorithm exhibits minor fluctuations (0.2-0.4) across the range of the assessed $\frac{\alpha_1}{\alpha_2}$ values without showing any particular trend. The lack of a clear correlation between RMSE and $\alpha_1/\alpha_2$ indicates that the algorithms tested operate with approximately the same accuracy independent of the assumed physical parameters. As shown in Figure 7, the GSBM is on average the most accurate in estimating the temperature over the domain. For $\alpha_1/\alpha_2 = 0.25$ the three algorithms display a similar RMSE. The observed oscillations in the RMSE values, especially in the GUM and GHYMLA algorithms, are mainly attributable to their proven sensitivity to the time increment. Note that for each algorithm the test has been conducted using the same $\Delta t$ for all the tested thermal diffusive gradients rather than tailoring the time increment to optimize accuracy for each specific $\frac{\alpha_1}{\alpha_2}$ value. By selecting the appropriate time step, the algorithms exhibit consistent accuracy regardless of the physical parameters of the domain.
A good accuracy in simulating the analytical solution is observed for all the assessed algorithms. The RMSE between $T_{\text{ref}}$ and $T_{\text{num}}$ is computed for an increasing number of particles $N_p$ as displayed in Figure 6. The GUM, GSBM, GHYMLA algorithms converge to the line $\text{RMSE} \propto N_p^{-0.5}$, i.e. the expected Monte Carlo convergence rate.

<figure><img src="image_8.png" /><figcaption>Figure 6: Case 2: RMSE between the analytical $T_{\text{ref}}$ and the numerical $T_{\text{num}}$ temperature obtained from GUM (red), GSBM (blue) and GHYMLA (black) for an increasing number of particles at time $t = 1$ s. The dashed line corresponds to $\text{RMSE} \propto N_p^{-0.5}$.</figcaption></figure>

Finally, we test the performance of the algorithms in estimating the temperature for different ratios between the thermal diffusivities of the two materials. The RMSE between the analytical and numerical temperature is computed for different ratio between the thermal diffusivity ($\alpha_1/\alpha_2$) in $R_1$ and $R_2$, as displayed in Figure 7. All algorithms exhibit RMSE fluctuating at a maximum of 10% with respect to the average without showing any particular trend. The lack of a clear correlation between RMSE and $\alpha_1/\alpha_2$ indicates that the algorithms tested operate with approximately the same accuracy independent of the assumed physical parameters. As shown in Figure 7, the GSBM is on average the most accurate in estimating the temperature over the domain. For $\alpha_1/\alpha_2 = 0.25$ the three algorithms display a similar RMSE.
<figure><img src="image_11.png" /><figcaption>Figure 7: Case 2: Comparison of the RMSE between $T_{\text{ref}}$ and $T_{\text{num}}$ derived from GUM, GSBM and GHYMLA for different $\alpha_1/\alpha_2$ values at time $t = 1$ s. The time step $\Delta t$ is assumed equal to $10^{-2}$ s for GSBM, $10^{-3}$ s for GUM and $10^{-4}$ s for GHYMLA.</figcaption></figure>

# 4 Conclusions

This work applies the modeling approach presented by [4] to simulate a thermal diffusion process in discontinuous media. The GUM, GHYMLA and GSBM algorithms are here employed to assess the conductive heat transport in one-dimensional bimaterial medium. The Uffink [30] and the HYMLA [12] methods have been developed to deal with bimaterial media, with a condition that the flux is continuous (Fick's law). For other interface conditions the Uffink and HYMLA methods cannot be applied as such. It is already noted by [21] that the SBM method could deal with a wide range of conditions. Here, we show through the GSBM method how to apply it to heat transport process. Although exact, the SBM method requires the simulation of complex distributions, which hinders its computational time. Fourier's law, which is a classical formulation of the evolution of the temperature, does not allow a direct application of these methods. We introduced a new formulation using a convenient transform and shown how to implement a Lagrangian method with a weight associated to each particle. The convenient transform is used to enable the application of the GUM, GSBM, and GHYMLA algorithms to heat transport problems. Our work provides a critical extension of available Lagrangian methodologies to deal with heat transfer processes in the presence of material interfaces. Based on the numerical results, the following major conclusions can be drawn:

* conservation of thermal energy is verified by the GUM, GSBM and GHYMLA algorithms which are all able to maintain steady-state thermal equilibrium condition;

* the time step magnitude strongly affects the accuracy of the algorithms in estimating the temperature over the domain in GUM and GHYMLA: the GUM demonstrates sensitivity to the time step only for $\Delta t > 0.1$ s, whereas GHYMLA exhibits a direct linear correlation (in log-log) between the root-mean-square error (RMSE) and the chosen $\Delta t$. The GSBM remains independent of $\Delta t$ over the span of three orders of magnitude;
<figure><img src="image_9.png" /><figcaption>Figure 7: Case 2: Comparison of the RMSE between $T_{\text{ref}}$ and $T_{\text{num}}$ derived from GUM (red circles), GSBM (blue squares) and GHYMLA (black diamonds) for different $\alpha_1/\alpha_2$ values at time $t = 1$ s. The time step $\Delta t$ is assumed equal to $10^{-2}$ s for GSBM, $10^{-3}$ s for GUM and $10^{-4}$ s for GHYMLA.</figcaption></figure>

# 4 Conclusions

This work applies the modeling approach presented in [2] to simulate a thermal diffusion process in discontinuous media. The GUM, GHYMLA and GSBM algorithms are here employed to assess the conductive heat transport in one-dimensional bimaterial medium. The Uffink [26] and the HYMLA [10] methods have been developed to deal with bimaterial media, with a condition that the flux is continuous (Fick's law). For other interface conditions the Uffink and HYMLA methods cannot be applied as such. It is already noted in [16] that the SBM method could deal with a wide range of conditions. Here, we show through the GSBM method how to apply it to heat transport process. Although exact, the SBM method requires the simulation of complex distributions, which hinders its computational time. Fourier's law, which is a classical formulation of the evolution of the temperature, does not allow a direct application of these methods. We introduced a new formulation using a convenient transform and shown how to implement a Lagrangian method with a weight associated to each particle.
* the tested algorithms demonstrate consistent accuracy regardless of the assumed physical parameters, such as thermal diffusivity contrast across the interface;

* the GUM algorithm strikes a favorable balance between accuracy and computational time, making it the preferred choice in the considered scenarios.

The key limitation of these algorithms is that they are valid in one-dimensional transport only. The extension of the method of images to multi-dimensional settings remains an open challenge, as it is inherently designed for one-dimensional transport problems. Consequently, the generalization of the proposed algorithms to higher dimensions is a promising avenue for future works. The extension of the proposed theory and algorithms to heat convection problems would be of a high interest for many applications and we plan to tackle this in a near future.

## Statements & Declarations

* **Funding:** 'The authors declare that no funds, grants, or other support were received during the preparation of this manuscript.'

* **Competing interests:** 'The authors have no relevant financial or non-financial interests to disclose.'

* **Authors' contributions:** 'All authors contributed equally to the study. All authors read and approved the final manuscript.'

* **Data availability:** 'The datasets generated during and/or analyzed during the current study are available from the corresponding author on reasonable request. The algorithms presented in this paper are available in a public repository https://gitlab.inria.fr/gsbm/gum-gsbm-ghymla'

## References

[1] Amin Amiri Delouei, Amin Emamian, Sajjad Karimnejad, Hasan Sajjadi, and Dengwei Jing. Two-dimensional temperature distribution in fgm sectors with the power-law variation in radial and circumferential directions. *Journal of Thermal Analysis and Calorimetry*, 144:611–621, 2021.

[2] Amin Amiri Delouei, Amin Emamian, Hasan Sajjadi, Meysam Atashafrooz, Yueming Li, Lian-Ping Wang, Dengwei Jing, and Gongnan Xie. A comprehensive review on multi-dimensional heat conduction of multi-layer and composite structures: Analytical solutions. *Journal of Thermal Science*, 30(6):1875 – 1907, 2021.

[3] Elisa Baioni, Antoine Lejay, Géraldine Pichot, and Giovanni Michele Porta. GUM-GSBM-GHYMLA. https://gitlab.inria.fr/gsbm/gum-gsbm-ghymla, 2024.

[4] Elisa Baioni, Antoine Lejay, Géraldine Pichot, and Giovanni Michele Porta. Modeling diffusion in discontinuous media under generalized interface conditions: theory and algorithms. *SIAM Journal on Scientific Computing*, 46(4):A2202–A2223, 2024.

[5] Horatio Scott Carslaw and John Conrad Jaeger. *Conduction of heat in solids*. Oxford University Press, 1959.
Our work provides a critical extension of available Lagrangian methodologies to deal with
heat transfer processes in the presence of material interfaces. Based on the numerical results,
the following major conclusions can be drawn:

* conservation of thermal energy is verified by the GUM, GSBM and GHYMLA algorithms which are all able to maintain steady-state thermal equilibrium condition;

* the time step magnitude strongly affects the accuracy of the algorithms in estimating the temperature over the domain in GUM and GHYMLA: the GUM demonstrates sensitivity to the time step only for $\Delta t > 0.1$ s, whereas GHYMLA exhibits a direct linear correlation (in log-log) between the root-mean-square error (RMSE) and the chosen $\Delta t$. The GSBM remains independent of $\Delta t$ over the span of three orders of magnitude;

* the tested algorithms demonstrate consistent accuracy regardless of the assumed physical parameters, such as thermal diffusivity contrast across the interface;

* the GUM algorithm strikes a favorable balance between accuracy and computational time, making it the preferred choice in the considered scenarios.

## References

[1] A. Amiri Delouei, A. Emamian, H. Sajjadi, M. Atashafrooz, Y. Li, L.-P. Wang, D. Jing, and G. Xie. A comprehensive review on multi-dimensional heat conduction of multi-layer and composite structures: Analytical solutions. *Journal of Thermal Science*, 30(6):1875–1907, 2021.

[2] E. Baioni, A. Lejay, G. Pichot, and G. M. Porta. Modeling diffusion in discontinuous media under generalized discontinuous interface conditions: theory and algorithms. working paper or preprint, July 2023.

[3] H. S. Carslaw and J. C. Jaeger. *Conduction of heat in solids*. Oxford University Press, 1959.

[4] W. Dahmen, T. Gotzen, S. Müller, and M. Rom. Numerical simulation of transpiration cooling through porous material. *International journal for numerical methods in fluids*, 76(6):331–365, 2014.

[5] F. de Monte. An analytic approach to the unsteady heat conduction processes in one-dimensional composite media. *International Journal of Heat and Mass Transfer*, 45(6):1333–1343, 2002.
[6] Wolfgang Dahmen, Thomas Gotzen, Siegfried Müller, and M Rom. Numerical simulation of transpiration cooling through porous material. *International journal for numerical methods in fluids*, 76(6):331–365, 2014.

[7] F. de Monte. An analytic approach to the unsteady heat conduction processes in one-dimensional composite media. *International Journal of Heat and Mass Transfer*, 45(6):1333–1343, 2002.

[8] Zhong-Shan Deng and Jing Liu. Monte Carlo method to solve multidimensional bioheat transfer problem. *Numerical Heat Transfer: Part B: Fundamentals*, 42(6):543–567, 2002.

[9] Muyassar E Ismaeel. Heat transfer in a square porous cavity with partial heating and cooling for opposite vertical walls. *Al-Rafidain Engineering Journal (AREJ)*, 19(5):107–121, 2011.

[10] J.B.J. Fourier. *The analytical theory of heat*. New York, Dover Publishers, 1955.

[11] Abdolhossein Haji-Sheikh and Ephraim M Sparrow. The floating random walk and its application to Monte Carlo solutions of heat equations. *SIAM Journal on Applied Mathematics*, 14(2):370–389, 1966.

[12] H. Hoteit, R. Mose, A. Younes, F. Lehmann, and Ph. 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] Arvydas J Janavičius and Sigita Turskiene. Modelling of nonlinear thermodiffusion for a spherically symmetric case. *East European Journal of Physics*, 1:13–19, 2021.

[14] Nima Fallah Jouybari, Majid Eshagh Nimvari, and Wennan Zhang. A comparative study of different heat transfer enhancement mechanisms in a partially porous pipe. *SN Applied Sciences*, 3:1–15, 2021.

[15] Zohaib Atiq Khan, Pablo Angel Garcia Salaberri, Thomas MM Heenan, Rhodri Jervis, Paul R Shearing, Dan Brett, Ali Elkamel, and Jeff T Gostick. Probing the structure-performance relationship of lithium-ion battery cathodes using pore-networks extracted from three-phase tomograms. *Journal of The Electrochemical Society*, 167(4):040528, 2020.

[16] Zohaib Atiq Khan, Tom Tranter, Mehrez Agnaou, Ali Elkamel, and Jeff Gostick. Dual network extraction algorithm to investigate multiple transport processes in porous materials: Image-based modeling of pore and grain scale processes. *Computers & Chemical Engineering*, 123:64–77, 2019.

[17] Timo Koch, Kilian Weishaupt, Johannes Müller, Bernhard Weigand, and Rainer Helmig. A (dual) network model for heat transfer in porous media: Toward efficient model concepts for coupled systems from fuel cells to heat exchangers. *Transport in Porous Media*, 140(1):107–141, 2021.

[18] Eric M. LaBolle, Graham E. Fogg, and Andrew F.B. 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.

[19] Eric M LaBolle, Jeremy Quastel, and Graham E Fogg. Diffusion theory for transport in porous media: Transition-probability densities of diffusion processes corresponding to advection-dispersion equations. *Water Resources Research*, 34(7):1685–1693, 1998.
[6] Z.-S. Deng and J. Liu. Monte Carlo method to solve multidimensional bioheat transfer problem. *Numerical Heat Transfer: Part B: Fundamentals*, 42(6):543–567, 2002.

[7] M. E Ismaeel. Heat transfer in a square porous cavity with partial heating and cooling for opposite vertical walls. *Al-Rafidain Engineering Journal (AREJ)*, 19(5):107–121, 2011.

[8] J. Fourier. *The analytical theory of heat*. New York, Dover Publishers, 1955.

[9] A. Haji-Sheikh and E. M. Sparrow. The floating random walk and its application to Monte Carlo solutions of heat equations. *SIAM Journal on Applied Mathematics*, 14(2):370–389, 1966.

[10] 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.

[11] A. J. Janavičius and S. Turskiene. Modelling of nonlinear thermodiffusion for a spherically symmetric case. *East European Journal of Physics*, 1:13–19, 2021.

[12] N. F. Jouybari, M. E. Nimvari, and W. Zhang. A comparative study of different heat transfer enhancement mechanisms in a partially porous pipe. *SN Applied Sciences*, 3:1–15, 2021.

[13] Z. A. Khan, P. A. G. Salaberri, T. M. Heenan, R. Jervis, P. R. Shearing, D. Brett, A. Elkamel, and J. T. Gostick. Probing the structure-performance relationship of lithium-ion battery cathodes using pore-networks extracted from three-phase tomograms. *Journal of The Electrochemical Society*, 167(4):040528, 2020.

[14] Z. A. Khan, T. Tranter, M. Agnaou, A. Elkamel, and J. Gostick. Dual network extraction algorithm to investigate multiple transport processes in porous materials: Image-based modeling of pore and grain scale processes. *Computers & Chemical Engineering*, 123:64–77, 2019.

[15] T. Koch, K. Weishaupt, J. Müller, B. Weigand, and R. Helmig. A (dual) network model for heat transfer in porous media: Toward efficient model concepts for coupled systems from fuel cells to heat exchangers. *Transport in Porous Media*, 140(1):107–141, 2021.

[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.
[20] A. Lejay and G. Pichot. Simulating diffusion processes in discontinuous media: Benchmark tests. *Journal of Computational Physics*, 314:384-413, 2016.

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

[22] Jiaqi Liu, Peixian Yu, Yingge Li, Chunhao Wan, and Dongxing Du. Numerical simulation on convective heat transfer characteristics in porous media based on the digital rock technology. *International Journal of Heat and Mass Transfer*, 196:123323, 2022.

[23] Ranit Monga, Rajdeep Deb, Daniel W Meyer, and Patrick Jenny. A probabilistic, flux-conservative particle-based framework for transport in fractured porous media. *Advances in Water Resources*, 172:104368, 2023.

[24] Libor Mrňa, Jan Řiháček, Martin Šarbort, and Petr Horník. Solar absorber with a structured surface—a way to increase efficiency. *Acta Polytechnica*, 59(2):134-143, April 2019.

[25] Benoit Noetinger, Delphine Roubinet, Anna Russian, Tanguy Le Borgne, Frederick Delay, Marco Dentz, Jean-Raynald De Dreuzy, and Philippe Gouze. Random walk methods for modeling hydrodynamic transport in porous and fractured media from pore to reservoir scale. *Transport in Porous Media*, 115:345-385, 2016.

[26] Matej Praprotnik, M Sterk, and Roman Trobec. Inhomogeneous heat-conduction problems solved by a new explicit finite difference scheme. *International Journal of Pure and Applied Mathematics*, 13:275-292, 2004.

[27] Morgan Sans, Olivier Farges, Vincent Schick, and Gilles Parent. Solving transient coupled conductive and radiative transfers in porous media with a Monte Carlo method: Characterization of thermal conductivity of foams using a numerical flash method. *International Journal of Thermal Sciences*, 179:107656, 2022.

[28] Ebrahim Shahraeeni and Dani Or. Thermo-evaporative fluxes from heterogeneous porous surfaces resolved by infrared thermography. *Water resources research*, 46(9), 2010.

[29] Thomas Sherman, Nicholas B. Engdahl, Giovanni Porta, and Diogo Bolster. A review of spatial Markov models for predicting pre-asymptotic and anomalous transport in porous and fractured media. *Journal of Contaminant Hydrology*, 236:103734, 2021.

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

[31] Nicole Vorhauer, Haashir Altaf, Evangelos Tsotsas, and Tanja Vidakovic-Koch. Pore network simulation of gas-liquid distribution in porous transport layers. *Processes*, 7(9):558, 2019.

[32] Kun Yang, Xin Li, Kangyuan Liu, and Jiabing Wang. Coupling effect of heat transfer in plate heat exchanger filled with porous media. *International Journal of Heat and Mass Transfer*, 182:121966, 2022.
[18] J. Liu, P. Yu, Y. Li, C. Wan, and D. Du. Numerical simulation on convective heat transfer characteristics in porous media based on the digital rock technology. *International Journal of Heat and Mass Transfer*, 196:123323, 2022.

[19] L. Mrňa, J. Řiháček, M. Šarbort, and P. Horník. Solar absorber with a structured surface—a way to increase efficiency. *Acta Polytechnica*, 59(2):134–143, Apr. 2019.

[20] B. Noetinger, D. Roubinet, A. Russian, T. Le Borgne, F. Delay, M. Dentz, J.-R. De Dreuzy, and P. Gouze. Random walk methods for modeling hydrodynamic transport in porous and fractured media from pore to reservoir scale. *Transport in Porous Media*, 115:345–385, 2016.

[21] M. Praprotnik, M. Sterk, and R. Trobec. Inhomogeneous heat-conduction problems solved by a new explicit finite difference scheme. *International Journal of Pure and Applied Mathematics*, 13:275–292, 2004.

[22] T. Sahabi and S. Balaska. Analytical one-dimensionless study of unsteady heat transfer in bilayer, and three-layer materials. 2020.

[23] M. Sans, O. Farges, V. Schick, and G. Parent. Solving transient coupled conductive and radiative transfers in porous media with a Monte Carlo method: Characterization of thermal conductivity of foams using a numerical flash method. *International Journal of Thermal Sciences*, 179:107656, 2022.

[24] E. Shahraeeni and D. Or. Thermo-evaporative fluxes from heterogeneous porous surfaces resolved by infrared thermography. *Water resources research*, 46(9), 2010.

[25] T. Sherman, N. B. Engdahl, G. Porta, and D. Bolster. A review of spatial markov models for predicting pre-asymptotic and anomalous transport in porous and fractured media. *Journal of Contaminant Hydrology*, 236:103734, 2021.

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

[27] Y. Xu, Y. Zheng, and J. Kou. Prediction of effective thermal conductivity of porous media with fractal-Monte Carlo simulations. *Fractals*, 22(03):1440004, 2014.

[28] K. Yang, X. Li, K. Liu, and J. Wang. Coupling effect of heat transfer in plate heat exchanger filled with porous media. *International Journal of Heat and Mass Transfer*, 182:121966, 2022.