# Dynamical analysis and numerical simulation of a reaction-diffusion model for microbial decomposition of organic matter in 3D soil structure

Mohammed Elghandouri, Klai Mouad, Khalil Ezzinbi, Olivier Monga

► To cite this version:

Mohammed Elghandouri, Klai Mouad, Khalil Ezzinbi, Olivier Monga. Dynamical analysis and nu-
merical simulation of a reaction-diffusion model for microbial decomposition of organic matter in 3D
soil structure. 2026. hal-05521937

HAL Id: hal-05521937

https://hal.science/hal-05521937v1

Preprint submitted on 21 Feb 2026

**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_1.png" /></figure>

Distributed under a Creative Commons CC BY 4.0 - Attribution - International License
Dynamical analysis and numerical simulation of a
reaction-diffusion model for microbial decomposition of
organic matter in 3D soil structure

Mohammed Elghandouri<sup>a,b,c,*</sup>, Klai Mouad<sup>b,c,d</sup>, Khalil Ezzinbi<sup>b,c</sup>, Olivier Monga<sup>c</sup>

<sup>a</sup>Inria, CNRS, Centrale Lyon, INSA Lyon, Université Claude Bernard Lyon 1, Université Jean Monnet, ICJ UMR5208, Villeurbanne, 69603, France

<sup>b</sup>LMDP Laboratory, Cadi Ayyad University, B.P: 2390, Marrakech, 40000, Morocco

<sup>c</sup>UMMISCO, IRD, Bondy, Paris, F-93143, France

<sup>d</sup>Centre for Remote Sensing Applications (CRSA), Mohammed VI Polytechnic University (UM6P), Lot 660, Hay Moulay Rachid, Ben Guerir, 43150, Morocco

Abstract

Microbial decomposition of organic matter in soil is a fundamental process in the global carbon cycle, directly influencing soil health, fertility, and greenhouse gas emissions. This paper presents a dynamic analysis and numerical simulation of a reaction-diffusion model that describes microbial decomposition of organic matter within a three dimensional soil structure. We investigate the interactions between Microbial Biomass (MB) and organic substrates, as well as the diffusion of various compounds through the soil matrix, using nonlinear parabolic partial differential equations. Our study provides proofs for the existence and uniqueness of solutions, as well as an analysis of asymptotic behavior. Notably, our investigation reveals the presence of a global attractor, where any solution, regardless of initial conditions, tends to converge. To illustrate the practical implications of our findings, we have developed a numerical tool to simulate the long-term behavior of the system with reasonable computational expense. This tool provides a visual proof of the global attractor for a validated set of biological parameters in a real sandy loam soil sample captured using 3D tomography imagery.

**Keywords:** Biological Dynamics in Soil, Non-linear Parabolic PDEs, Global Attractor, Pore Space Modelling, Simulation in Complex Shapes.

2020 MSC: 35A01, 35B41, 37C70, 65L12.

<footnote>*Corresponding author</footnote>

<footnote>Email addresses: m.elghandouri.ced@uca.ac.ma (Mohammed Elghandouri), mouad.klai-ext@um6p.ma (Klai Mouad), ezzinbi@uca.ac.ma (Khalil Ezzinbi), olivier.monga@ird.fr (Olivier Monga)</footnote>
Introduction

Microbial activity in soil is vital for preserving soil health and ensuring ecosystem functionality. Microorganisms play key roles in nutrient cycling, breaking down organic matter, and fostering plant growth [1, 2]. The vast diversity and intricate interactions within soil present a challenge to fully grasp its complexity. Furthermore, soil's heterogeneous and multifaceted nature makes directly observing and measuring microbial activity difficult [3]. Traditional methods, such as culturing and microscopy, are not only time-consuming but also offer limited insights. Hence, there is a pressing need for alternative methods for studying soil microbial activity [4]. Over the past two decades, numerous studies have concentrated on modeling the dynamics within the complex soil structure. These efforts seek to mimic how soil components move and interact in porous media and fractured environments. Several significant contributions have emerged from these studies, with models typically integrating transport mechanisms with sophisticated reaction (transformation) processes to simulate the activity of microorganisms and compounds within the soil matrix. Significant attention has been given to computational models, focusing on the development of numerical tools to simulate transformation and diffusion processes in the complex pore space of soil. Various methodologies have been employed to describe chemical transport phenomena, primarily the Lattice-Boltzmann Method (LBM) [7] and pore network geometrical models (also known as morphological models) [9, 10, 12, 13]. Another work, presented in [11], proposes a region-based representation in which the pore space is partitioned into connected subsets of voxels derived from the curve skeleton, and an attributed adjacency graph is constructed to model pore connectivity and transport phenomena. On the other hand, the mathematical analysis of reaction-diffusion models has only been investigated in a limited number of studies. One notable example is found in [14], where the authors used nonlinear parabolic Partial Differential Equations (PDEs) to model the microbial decomposition of organic matter in soil. This approach allowed for a detailed representation of the biochemical processes occurring within the soil matrix. To solve these complex equations, the authors employed FreeFEM++, a powerful finite element software, to simulate the model based on the variational formulation of the PDEs. However, further research is needed to analyse and refine these mathematical models.

In the present paper, we investigate a nonlinear parabolic PDE model describing
microbial activity in soil and provide a formal analysis of the model, proving the
existence and uniqueness of a global solution. Then, we prove the existence of
a global attractor, where any solution starting from any initial distribution tends
to converge. Furthermore, we provide a numerical tool to simulate the long-term
behavior of the system in a reasonable time using an optimal sphere network that
approximates the pore space of soil. Our research begins by describing a biological
model that represents the dynamic processes occurring within the soil. This model
encapsulates key biological parameters and system behaviors, shedding light on
organic matter decomposition by microorganisms and the production of mineralized
CO₂, and taking into account the diffusion of different compounds through the soil matrix.

To rigorously analyze this biological system, we employ PDEs as our mathematical tool of choice. Specifically, we formulate a set of parabolic nonlinear PDEs that elegantly capture the diffusion and transformation processes expressing the degradation of organic matter by microorganisms taking place within the three-dimensional (3D) soil structure. These equations are accompanied by Neumann boundary conditions, reflecting the system's behavior at its boundaries. Before diving into the core of our analysis, we provide an overview of the mathematical tools and concepts essential to our proof of the global attractor for the PDE model [17, 22, 24, 18, 23]. These preliminary insights based on dynamical systems theory set the stage for the subsequent sections, where we unveil the intriguing properties of the system. One of our primary objectives is to establish the global existence of solutions to the PDE model. We demonstrate that the problem possesses a unique positive mild solution defined over an unbounded time interval, $[0, +\infty[$. This foundational result forms the basis for our exploration of the system's long-term behavior. At the core of our research, we reveal the presence of a global attractor which is a compact, connected and invariant subset of the positive points comprising the system's phase space. Importantly, this global attractor possesses the remarkable capacity to attract and encapsulate all bounded subsets of the phase space, thereby affording invaluable insights into the long-term behavior of the system. To enhance the robustness of our theoretical findings, we employ a geometric approach to model the intricate pore space within the soil structure. This methodology involves adapting the PDE model to fit the novel framework of the pore network model, a technique that has been previously validated [10, 12, 13]. The approach exhibits exceptional capabilities in simulating the long-term behavior of the system and is considered as a powerful alternative to the LBM in simulating diffusion in complex geometries. In our study, we employ this method to explore a range of realistic parameters and scenarios using a real soil sample obtained through computed tomographic imaging.

The description of this work is organized as follows. In Section 1, we present our biological model. In Section 2, we introduce our model using PDEs. Section 3 revisits essential properties related to global attractors for partial differential equations and discusses properties crucial for proving our main results. In Section 4, we delve into the global existence of our PDEs system. Section 5 establishes the existence of a global attractor for our model. Section 6 is dedicated to presenting numerical simulations for visualizing the attractor. Finally, we showcase and discuss our results in Section 7.

# 1. General description of the biological model

Our objective is to simulate and model the process of organic carbon mineralization by microorganisms within the soil pore space. We improve the methodology outlined in [9, 19, 14] to simulate the organic matter decomposition by microorganisms in the soil pore space by solving a PDEs system. Microorganisms degrade
organic matter in accordance with the laws of supply and demand. The method presented in [14] is based on the mass conservation law to obtain a reaction-diffusion model that describes the diffusion of organic matter. In our approach, we assume that the decomposition of soil organic matter involves five key biological components.

Microbial activity is governed by the diffusion of various compounds within the soil's pores and by transformation processes. We model the dynamics of microbial decomposition using five key compounds, as illustrated in figure 1:

*   **Microbial Biomass (MB):** Represents the mass of microorganisms in the sample.

*   **Carbon Dioxide ($CO_2$):** Indicates the carbon dioxide produced by the breathing of the micro-organisms. The quantity (mass) of $CO_2$ emitted per time unit is proportional to the biomass and is defined by the respiration rate multiplied by the biomass.

*   **Fresh Organic Matter (FOM):** Derived from recently added or deposited plant and easily decomposable.

*   **Soil Organic Matter (SOM):** Consists of various organic compounds in different stages of decomposition, originating from biomass turnover and less accessible to decomposition.

*   **Dissolved Organic Matter (DOM):** Refers to organic compounds dissolved in soil water, originating from hydrolysis of FOM and SOM and biomass recycling, available for microbial uptake or transport within the pore space.

The DOM arises from the decomposition of both slow-decomposing SOM and fast-decomposing FOM. Microorganisms grow through the assimilation of DOM and respire by producing $CO_2$. Upon death, they are transformed back into organic matter. The biological process is summarized in Figure 1.

## 2. Mathematical formulation of the model

Let $\Omega \subset \mathbb{R}^3$ be the domain representing the 3D soil pore space, $t \ge 0$ the time given, and $x = (x_1, x_2, x_3)^\top \in \Omega$ be a point of the pore space. We denote by $b(x, t)$, $n(x, t)$, $m_1(x, t)$, $m_2(x, t)$, and $c(x, t)$ the densities of the MB, the DOM, the SOM, the FOM, and the $CO_2$ at position $x \in \Omega$ and time $t \ge 0$ respectively.

Let $\mathcal{V}$ be an individual volume in $\Omega$. The evolution of the density of the MB in $\mathcal{V}$ depends on the diffusion of the microorganisms, the growth of the microorganisms, the mortality of the microorganisms, and the breathing of the microorganisms. It is assumed that the MB consumes the DOM to grow. The variation of the MB density is described by the following equation:

$$ \frac{\partial b}{\partial t} = D_b \Delta b + \frac{Kn}{K_b + n} b - (\mu + \eta)b, $$
<figure><img src="image_2.png" /><figcaption>Figure 1: Biological processes simulated inside one spatial unit and the interactions with external units.</figcaption></figure>

where $D_b$ is the microbial decomposer's diffusion coefficient, $K$ is the maximal growth rate, $K_b$ is the half-saturation constant, $\mu$ is the mortality rate, and $\eta$ is the breathing rate.

The evolution of the density of the DOM in $V$ depends on the DOM diffusion, the assimilation of the DOM by the MB decomposer's, the MB decomposer's mortality, and the transformation of the SOM and the FOM. The equation that describes the DOM density variation is the following:

$$ \frac{\partial n}{\partial t} = D_n \Delta n - \frac{Kn}{K_b + n} b + c_1 m_1 + c_2 m_2 + \rho \mu b, $$

where $D_n$ is the DOM diffusion coefficient, $c_1$ is the maximal transformation rate of the SOM, $c_2$ is the maximal transformation rate of the FOM, and $\rho\mu$ is the transformation rate of the deceased MB decomposer's into the DOM.

The variation in SOM quantity arises from the conversion of a portion of the SOM into the DOM and from the mortality of MB. Thus, the evolution of the density of the SOM in $V$ is presented by the following equation:

$$ \frac{\partial m_1}{\partial t} = -c_1 m_1 + (1 - \rho) \mu b. $$

The evolution of the density of the FOM in $V$ is presented by the following equation:

$$ \frac{\partial m_2}{\partial t} = -c_2 m_2. $$

The evolution of the CO₂ density in $V$ depends on its diffusion and its production by the MB decomposer's during breathing. The CO₂ quantity variation is represented by the following equation:

$$ \frac{\partial c}{\partial t} = D_c \Delta c + \eta b, $$
where $D_c$ is the CO$_2$ diffusion coefficient.

Therefore, the microbial decomposition of organic matter in soil is described by the following PDEs system:

$$ 
\left\{
\begin{aligned}
\frac{\partial b}{\partial t} &= D_b \Delta b + \frac{Kn}{K_b + n} b - (\mu + \eta)b, \\
\frac{\partial n}{\partial t} &= D_n \Delta n - \frac{Kn}{K_b + n} b + c_1 m_1 + c_2 m_2 + \rho \mu b, \\
\frac{\partial m_1}{\partial t} &= -c_1 m_1 + (1 - \rho) \mu b, \\
\frac{\partial m_2}{\partial t} &= -c_2 m_2, \\
\frac{\partial c}{\partial t} &= D_c \Delta c + \eta b.
\end{aligned}
\right.
\quad (1) $$

On the border $\partial\Omega$ of $\Omega$, we use Neumann boundary conditions, which means that the flow is null on $\partial\Omega$ for all variables, i.e,

$$ \frac{\partial b}{\partial \vec{v}} = \frac{\partial n}{\partial \vec{v}} = \frac{\partial m_1}{\partial \vec{v}} = \frac{\partial m_2}{\partial \vec{v}} = \frac{\partial c}{\partial \vec{v}} = 0, \quad \text{on } \partial\Omega \times \mathbb{R}^+, \quad (2) $$

where $\vec{v}$ denotes the outward unit normal vector.

For the initial conditions in $\Omega$, we denote by

$$ (b(0, \cdot), n(0, \cdot), m_1(0, \cdot), m_2(0, \cdot), c(0, \cdot)) = (b_0, n_0, m_{1,0}, m_{2,0}, c_0). \quad (3) $$

To simplify (1)-(3) into a vector form equation, we define: the state system $u$ as:

$$ u = (u_1, u_2, u_3, u_4, u_5)^{\top} = (b, n, m_1, m_2, c)^{\top}. $$

The initial state $u_0$ as:

$$ u_0 = (b_0, n_0, m_{1,0}, m_{2,0}, c_0)^{\top}. $$

The diffusion coefficients matrix $\underline{D}$ as:

$$ \underline{D} = \begin{pmatrix} D_b & 0 & 0 & 0 & 0 \\ 0 & D_n & 0 & 0 & 0 \\ 0 & 0 & 0 & 0 & 0 \\ 0 & 0 & 0 & 0 & 0 \\ 0 & 0 & 0 & 0 & D_c \end{pmatrix}. $$

The reaction diffusion terms vector $F(u)$ is given by:

$$ F(u) = (F_1(u), F_2(u), F_3(u), F_4(u), F_5(u))^\top, $$
where

$$F_1(u) = \frac{K u_2}{K_b + u_2} u_1 - (\mu + \eta) u_1,$$

$$F_2(u) = -\frac{K u_2}{K_b + u_2} u_1 + c_1 u_3 + c_2 u_4 + \rho \mu u_1,$$

$$F_3(u) = -c_1 u_3 + (1 - \rho) \mu u_1,$$

$$F_4(u) = -c_2 u_4,$$

$$F_5(u) = \eta u_1.$$

Finally, equation (1)-(3) can be reduced to:

$$\left\{
\begin{array}{ll}
\displaystyle \frac{\partial u}{\partial t} = \underline{D} \Delta u + F(u), & \text{in } \Omega \times \mathbb{R}^{+}, \\
\displaystyle \frac{\partial u}{\partial \vec{\nu}} = 0, & \text{on } \partial \Omega \times \mathbb{R}^{+}, \\
u(0) = u_0, & \text{in } \Omega.
\end{array}
\right.
\qquad (4)$$

## 2.1. Why do we need to look for an attractor?

Solving the equation $F(u) = 0$ allows us to determine the set of constant equilibrium points of equation (4), which is given by:

$$S = \{(0, u_2^*, 0, 0, u_5^*) : u_2^*, u_5^* \in \mathbb{R}^+\}.$$

This notation indicates that the equilibrium points are characterized by specific values for the variables $u_2^*$ and $u_5^*$, both of which belong to the set of non-negative real numbers. An important observation is that these equilibrium points are not isolated. This lack of isolation has implications for our analytical methods. Specifically, it means that we cannot employ techniques like linearization or Lyapunov's method, which are often used to analyze the behavior of systems near isolated equilibrium points. Furthermore, we do not expect to observe convergence to the equilibrium points, which makes understanding the system's behavior as it approaches equilibrium challenging. Given these challenges, our emphasis lies in establishing the existence of a global attractor for equation (4).

## 3. Basic working tools

In the sequel, we recall some useful properties on global attractor for partial differential equations. For more details, we refer to [17].

**Definition 3.1.** [6, Definition 2.4, page 143] A semiflow $(\mathcal{U}(t))_{t \ge 0}$ on a complete metric space $\mathcal{X}$ is a one-parameter family of maps $\mathcal{U}(t) : \mathcal{X} \to \mathcal{X}$, parameter $t \in \mathbb{R}^{+}$, such that
i) $\mathcal{U}(0) = \mathrm{id}_\mathcal{X}$ (id$_\mathcal{X}$ the identity map on $\mathcal{X}$).

ii) $\mathcal{U}(t+s) = \mathcal{U}(t)\mathcal{U}(s)$ for $t, s \ge 0$.

iii) $t \to \mathcal{U}(t)x$ is continuous for all $x \in \mathcal{X}$.

**Definition 3.2.** [17, Section 3.1, page 9] Let $(\mathcal{U}(t))_{t \ge 0}$ be a semiflow on a complete metric space $\mathcal{X}$.

i) We say that $B \subset \mathcal{X}$ is invariant under $(\mathcal{U}(t))_{t \ge 0}$ if $\mathcal{U}(t)B = B$ for all $t \ge 0$.

ii) We say that $(\mathcal{U}(t))_{t \ge 0}$ is bounded, if it takes bounded sets into bounded sets.

**Definition 3.3.** [17, Section 2.4, page 16] Let $(\mathcal{U}(t))_{t \ge 0}$ be a semiflow on a complete metric space $\mathcal{X}$. We say that $(\mathcal{U}(t))_{t \ge 0}$ is point dissipative if there is a bounded set $B \subset \mathcal{X}$ which attracts each point of $\mathcal{X}$ under $(\mathcal{U}(t))_{t \ge 0}$.

**Theorem 3.1.** [17] Let $(\mathcal{U}(t))_{t \ge 0}$ be a semiflow on a complete space $\mathcal{X}$ equipped with the metric $\| \cdot \|$. Assume that there exists $L > 0$ such that for all $x_0 \in \mathcal{X}$, $\limsup_{t \to +\infty} \|\mathcal{U}(t)x_0\| \le L$, then $(\mathcal{U}(t))_{t \ge 0}$ is point dissipative.

**Definition 3.4.** [17, Section 3.1, page 36] Let $(\mathcal{U}(t))_{t \ge 0}$ be a semiflow on a complete metric space $(\mathcal{X}, d)$; $d$ is the metric in $\mathcal{X}$. A subset $B$ of $\mathcal{X}$ is said to attract $C \subset \mathcal{X}$ if $d(\mathcal{U}(t)C, B) \to 0$ as $t \to +\infty$, where

$$d(\mathcal{U}(t)C, B) := \inf\{d(\mathcal{U}(t)x, y); x \in B \text{ and } y \in C\}.$$

**Definition 3.5.** [17, Section 2.4, page 39] Let $\mathcal{X}$ be a complete metric space and $(\mathcal{U}(t))_{t \ge 0}$ be a semiflow on $\mathcal{X}$. Let $\mathcal{A}$ be a subset of $\mathcal{X}$. Then, $\mathcal{A}$ is called a global attractor for $(\mathcal{U}(t))_{t \ge 0}$ in $\mathcal{X}$ if $\mathcal{A}$ is a compact, connected, and invariant set of $\mathcal{X}$ that attracts every bounded set of $\mathcal{X}$.

The following Theorem shows the existence of global attractor for point dissipative semiflows.

**Theorem 3.2.** [17, Theorem 4.1.2, page 63] Let $(\mathcal{U}(t))_{t \ge 0}$ be a semiflow on a complete metric space $\mathcal{X}$. Suppose that $\mathcal{U}(t)$ is bounded for each $t \ge 0$. If $(\mathcal{U}(t))_{t \ge 0}$ is point dissipative, then it has a connected global attractor.

In the sequel, we need the following Lemma.

**Lemma 3.1.** [15, Theorem 11.3 with $h = 0$, page 99] Let $v$ be a real, continuous and nonnegative function such that

$$v(t) \le c + \int_{t_0}^{t} w(t, s)v(s)ds, \quad \text{for } t \ge t_0,$$

where $c > 0$, $w(t, s)$ is continuously differentiable in $t$ and continuous in $s$ with $w(t, s) \ge 0$ for $t \ge s \ge t_0$. Then,

$$v(t) \le c \exp \left( \int_{t_0}^{t} \left[ w(s, s) + \int_{t_0}^{s} \frac{\partial w(s, r)}{\partial s} dr \right] ds \right), \quad \text{for } t \ge t_0.$$
### 4. Global existence for equation (4)

The space $(\mathbb{L}^2(\Omega))^5$ endowed with the norm $\|\cdot\|_{(\mathbb{L}^2(\Omega))^5} = \sum \|\cdot\|_{\mathbb{L}^2(\Omega)}$ issue from the inner product $(\cdot, \cdot)$ defined by

$$ (v, w) = \sum_{i=1}^{5} \langle v_i, w_i \rangle_{\mathbb{L}^2(\Omega)} = \sum_{i=1}^{5} \int_{\Omega} v_i(x) w_i(x) dx, $$

for $v = (v_1, v_2, v_3, v_4, v_5)$, $w = (w_1, w_2, w_3, w_4, w_5) \in (\mathbb{L}^2(\Omega))^5$ is a Hilbert space.

In order to rewrite equation (4) in an abstract form, we define the operator $(A, D(A))$ by

$$ \left\{ \begin{array}{l} D(A) = (H_0^2(\Omega))^2 \times (\mathbb{L}^2(\Omega))^2 \times H_0^2(\Omega), \\ Av = \underline{D}\Delta v, \text{ for } v \in D(A), \end{array} \right. $$

where

$$ H_0^2(\Omega) = \left\{ w \in \mathbb{L}^2(\Omega) : \text{ such that } \nabla w, \Delta w \in \mathbb{L}^2(\Omega), \text{ and } \frac{\partial w}{\partial \overline{v}} = 0 \right\}. $$

Then, equation (4) can be transformed to the following reaction-diffusion form:

$$ \left\{ \begin{array}{l} u'(t) = Au(t) + F(u(t)), \text{ for } t \ge 0, \\ u(0) = u_0. \end{array} \right. \qquad (5) $$

Let $\lambda > 0$ be fixed. Equation (5) can be written in the following equivalent form:

$$ \left\{ \begin{array}{l} u'(t) = (A - \lambda I)u(t) + G(u(t)), \text{ for } t \ge 0, \\ u(0) = u_0, \end{array} \right. \qquad (6) $$

where $G(v) = F(v) + \lambda v$ for $v \in (\mathbb{L}^2(\Omega))^5$.

In the next, we denote by

$$ X := (\mathbb{L}^2(\Omega))^5 \text{ endowed with the norm } \|\cdot\| := \|\cdot\|_{(\mathbb{L}^2(\Omega))^5}, $$

and

$$ \Lambda_+ := \{ v = (v_1, v_2, v_3, v_4, v_5)^\top \in X : v_i \ge 0, \text{ for } i = 1, \dots, 5 \}. $$

**Remark 4.1.** We have

$$ G(v) = F(v) + \lambda v = \begin{pmatrix} \frac{K v_2}{K_b + v_2} v_1 + (\lambda - (\mu + \eta)) v_1 \\ c_1 v_3 + c_2 v_4 + \left(\lambda - \frac{K v_1}{K_b + v_2}\right) v_2 + \rho \mu v_1 \\ (\lambda - c_1) v_3 + (1 - \rho) \mu v_1 \\ (\lambda - c_2) v_4 \\ \eta v_1 + \lambda v_5 \end{pmatrix}, $$

for each $v \in \Lambda_+$. By Definition 4.1 in [8], we can show that $G$ is essentially nonnegative on $\Lambda_+$ for a sufficiently large choice of $\lambda > 0$.
**Lemma 4.1.** There exists a positive constant $C_1 > 0$ such that $||F(v)|| \le C_1||v||$ for each $v \in \Lambda_+$.

**Proof:** Let $v = (v_1, v_2, v_3, v_4, v_5)^\top \in \Lambda_+$. Since,

$$|F_1(v)(x)|^2 \le 2(K^2 + (\mu + \eta)^2) v_1(x)^2;$$

$$|F_2(v)(x)|^2 \le 4[(K^2 + (\rho\mu)^2)v_1(x)^2 + c_1^2v_3(x)^2 + c_2^2v_4(x)^2];$$

$$|F_3(v)(x)|^2 \le 2[c_1^2v_3(x)^2 + ((1 - \rho)\mu)^2v_1(x)^2];$$

$$|F_4(v)(x)|^2 \le c_2^2v_4(x)^2;$$

$$|F_5(v)(x)|^2 \le \eta^2v_1(x)^2.$$

It follows that,

$$\begin{aligned}
\|F(v)\| &= \sum_{i=1}^{5} \|F_i(v)\|_{L^2(\Omega)} \\
&= \sum_{i=1}^{5} \left( \int_{\Omega} |F_i(v)(x)|^2 dx \right)^{1/2} \\
&= \left( \int_{\Omega} 2(K^2 + (\mu + \eta)^2) v_1(x)^2 dx \right)^{1/2} \\
&\quad + \left( \int_{\Omega} 4[(K^2 + (\rho\mu)^2)v_1(x)^2 + c_1^2v_3(x)^2 + c_2^2v_4(x)^2] dx \right)^{1/2} \\
&\quad + \left( \int_{\Omega} 2[c_1^2v_3(x)^2 + ((1-\rho)\mu)^2v_1(x)^2] dx \right)^{1/2} + \left( c_2^2 \int_{\Omega} v_4(x)^2 dx \right)^{1/2} \\
&\quad + \left( \eta^2 \int_{\Omega} v_1(x)^2 dx \right)^{1/2} \\
&= \sqrt{2}(K^2 + (\mu + \eta)^2)^{1/2} \|v_1\|_{L^2(\Omega)} \\
&\quad + 2(K^2\|v_1\|_{L^2(\Omega)}^2 + c_1^2\|v_3\|_{L^2(\Omega)}^2 + c_2^2\|v_4\|_{L^2(\Omega)}^2 + (\rho\mu)^2\|v_1\|_{L^2(\Omega)}^2)^{1/2} \\
&\quad + \sqrt{2}(c_1^2\|v_3\|_{L^2(\Omega)}^2 + ((1-\rho)\mu)^2\|v_1\|_{L^2(\Omega)}^2)^{1/2} + c_2\|v_4\|_{L^2(\Omega)} + \eta\|v_1\|_{L^2(\Omega)}.
\end{aligned}$$

We use the fact that $\sqrt{a+b} \le \sqrt{a} + \sqrt{b}$ for $a, b \ge 0$, we get that

$$\begin{aligned}
\|F(v)\| &\le \left( \sqrt{2(K^2 + (\mu + \eta)^2)} + \eta + 2K + \sqrt{2}((1 - \rho)\mu)^2 \right) \|v_1\|_{L^2(\Omega)} \\
&\quad + \sqrt{2}(\sqrt{2} + 1) c_1 \|v_3\|_{L^2(\Omega)} + 3c_2 \|v_4\|_{L^2(\Omega)} \\
&\le C_1 \|v\|,
\end{aligned}$$
where

$$C_1 = \max \left[ 1, \sqrt{2(K^2 + (\mu + \eta)^2)} + \eta + 2K + \sqrt{2((1-\rho)\mu)^2}, \sqrt{2}(\sqrt{2}+1)c_1, 3c_2 \right]. \quad \square$$

**Lemma 4.2.** For each $r > 0$, there exists a positive constant $C_2(r) > 0$ such that

$$\|F(u) - F(v)\| \le C_2(r)\|u - v\|, \text{ for } u, v \in B(0, r) \cap \Lambda_+,$$

where $B(0, r)$ is the open ball (in $X$) of center 0 and radius $r$.

**Proof:** The proof follows the fact that $F$ is $C^1$ from $\Lambda^+$ to $X$. $\square$

**Theorem 4.1.** Operator $(A, D(A))$ generates a $C_0$-semigroup $(T(t))_{t \ge 0}$ of contractions on $X$.

**Proof:** We show that $-A$ is maximal monotone. Let $D_1 = D_b$, $D_2 = D_n$, $D_3 = D_4 = 0$, and $D_5 = D_c$. Let $v = (v_1, \dots, v_5)^\top \in D(A)$, then

$$\begin{align*}
\langle v, -Av \rangle &= \sum_{i=1}^{5} \int_{\Omega} -D_i \Delta v_i(x) v_i(x) dx \\
&= \sum_{i=1}^{5} D_i \int_{\Omega} |\nabla v_i(x)|^2 dx \\
&= \sum_{i=1}^{5} D_i \|\nabla v_i\|_{L^2(\Omega)}^2.
\end{align*}$$

Hence, $-A$ is monotone. Since $rg(I - D_i\Delta) = L^2(\Omega)$ for each $i = 1, 2, \dots, 5$ (where $rg(I - D_i\Delta)$ denotes the range of the operator $I - D_i\Delta$), it follows that for every $v = (v_1, v_2, v_3, v_4, v_5)^\top \in X$, there exists $w = (w_1, w_2, w_3, w_4, w_5)^\top \in D(A)$ such that $(I - A)w = v$. Thus, $rg(I - A) = X$, that is $-A$ is maximal. Consequently, $(A, D(A))$ is an infinitesimal generator of a $C_0$-semigroup $(T(t))_{t \ge 0}$ of contractions on $X$ [21]. $\square$

**Theorem 4.2.** [16] The $C_0$-semigroup $(T(t))_{t \ge 0}$, is positive, namely $T(t)\Lambda_+ \subseteq \Lambda_+$ for all $t \ge 0$.

**Remark 4.2.** For each $\lambda > 0$, operator $(A - \lambda I, D(A))$ generates a $C_0$-semigroup $(S(t))_{t \ge 0}$ given by

$$S(t) = e^{-\lambda t} T(t), \text{ for } t \ge 0.$$

**Definition 4.1.** A continuous function $u : \mathbb{R}^+ \to X$ is said to be a mild solution of equation (6) if
$$u(t) = S(t)u_0 + \int_0^t S(t-s)G(u(s))ds, \text{ for } t \ge 0.$$

We recall the following theorem. The fact that $f$ is at most affine ensures global existence.

**Theorem 4.3.** [21, Theorem 1.4, page 185] Let $f: X \to X$ be locally Lipschitz continuous and at most affine. If $B$ is the infinitesimal generator of a $C_0$-semigroup $(T(t))_{t \ge 0}$ on $X$, then for every $v_0 \in X$, the initial value problem

$$\begin{cases} v'(t) = Bv(t) + f(v(t)), & t \ge 0, \\ v(0) = v_0, \end{cases}$$

has a unique mild solution $v$ on $[0, +\infty[$.

The following Theorem is the main result in this section.

**Theorem 4.4.** Assume that $u_0 \in \Lambda_+$. Equation (6) has a unique positive mild solution defined on $[0, +\infty[$.

**Proof:** The existence and uniqueness of the solution over the interval $[0, +\infty[$ can be readily established based on the findings of Lemma 4.1, Lemma 4.2, and Theorem 4.3. Additionally, as noted in Remark 4.1, we can choose a sufficiently large positive value for $\lambda$ to ensure that function $G$ is essentially nonnegative on $\Lambda_+$. Subsequently, proving the positivity of the solution reduces to Proposition 4.1 in [8] and Theorem 4.2. $\square$

**5. Global attractor for equation (4)**

Let us define the family of mappings $\{\mathcal{U}(t)\}_{t \ge 0}$ from $\Lambda_+$ to $\Lambda_+$ by

$$\mathcal{U}(t)u_0 := u(t, u_0) := u(t), \quad \text{for all } t \ge 0,$$

where $u(t)$ denotes the mild solution of equation (6) at time $t$ corresponding to the initial condition $u_0$. It follows immediately that $\mathcal{U}(0)u_0 = u_0$ and $t \to \mathcal{U}(t)u_0$ is continuous for all $u_0 \in \Lambda_+$. Since (6) is autonomous, the property

$$\mathcal{U}(t+s) = \mathcal{U}(t)\mathcal{U}(s), \quad \text{for } t \ge s \ge 0,$$

is a consequence of the uniqueness of the solution. As a consequence, $(\mathcal{U}(t))_{t \ge 0}$ is a semiflow on $\Lambda_+$. The following result shows that $\mathcal{U}(t)$ is a bounded map for each $t \ge 0$.

**Proposition 5.1.** There exists $\tilde{C}_1 > 1$ such that $\|\mathcal{U}(t)u_0\| \le \tilde{C}_1\|u_0\|$ for $t \ge 0$ and all $u_0 \in \Lambda_+$.
**Proof:** Let $u_0 \in \Lambda_+$ be fixed. Using Lemma 4.1, we show that for each $t \ge 0$, we have

$$ \begin{aligned} \|\mathcal{U}(t)u_0\| &\le e^{-\lambda t}\|u_0\| + \int_0^t e^{-\lambda(t-s)}\|G(\mathcal{U}(s)u_0)\|ds \\ &\le e^{-\lambda t}\|u_0\| + \int_0^t e^{-\lambda(t-s)}(\lambda + C_1)\|\mathcal{U}(s)u_0\|ds \\ &\le \|u_0\| + \int_0^t w(t,s)\|\mathcal{U}(s)u_0\|ds, \end{aligned} $$

where

$$ w(t,s) = e^{-\lambda(t-s)}(\lambda + C_1), \quad \text{for } t \ge s \ge 0. $$

By Lemma 3.1, we obtain that

$$ \begin{aligned} \|\mathcal{U}(t)u_0\| &\le \exp \left( \int_0^t \left[ (\lambda + C_1) + \int_0^s \frac{\partial w(s,r)}{\partial s} dr \right] ds \right) \|u_0\| \\ &= \exp \left( \int_0^t (\lambda + C_1) \left[ 1 - \lambda \int_0^s e^{-\lambda(s-r)} dr \right] ds \right) \|u_0\| \\ &= \exp \left( \int_0^t (\lambda + C_1) \left[ 1 - (1 - e^{-\lambda s}) \right] ds \right) \|u_0\| \\ &= \exp \left( \int_0^t (\lambda + C_1) e^{-\lambda s} ds \right) \|u_0\| \\ &= \exp \left( \frac{(\lambda + C_1)}{\lambda} (1 - e^{-\lambda t}) \right) \|u_0\| \\ &\le \tilde{C}_1 \|u_0\|, \end{aligned} $$

where

$$ \tilde{C}_1 = \exp \left( 1 + \frac{C_1}{\lambda} \right). $$

□

**Theorem 5.1.** The semiflow $(\mathcal{U}(t))_{t \ge 0}$ is point dissipative, i.e, there exists $\tilde{C}_2 > 0$ such that $\limsup_{t \to +\infty} \|\mathcal{U}(t)u_0\| \le \tilde{C}_2$ for all $u_0 \in \Lambda_+$.

**Proof:** Let $u_0 \in \Lambda_+$ be fixed, and $t_0 := t_0(u_0) > 0$ be such that

$$ \Theta(t) := e^{-\lambda t} \left[ 1 + \frac{(\lambda + C_1)\tilde{C}_1}{\lambda} (e^{\lambda t_0} - 1) \right] \|u_0\| \le 1, \quad \text{for } t > t_0. $$
Then, for each $t > t_0$, we have

$$ \|\mathcal{U}(t)u_0\| \le e^{-\lambda t} \|u_0\| + \int_0^{t_0} e^{-\lambda(t-s)}(\lambda + C_1) \|\mathcal{U}(s)u_0\| ds \\ + \int_{t_0}^t e^{-\lambda(t-s)}(\lambda + C_1) \|\mathcal{U}(s)u_0\| ds. $$

By Proposition 5.1, we obtain that

$$ \begin{aligned} \|\mathcal{U}(t)u_0\| &\le e^{-\lambda t} \|u_0\| + (\lambda + C_1) \tilde{C}_1 \left( \frac{e^{-\lambda(t-t_0)} - e^{-\lambda t}}{\lambda} \right) \|u_0\| \\ &\quad + \int_{t_0}^t e^{-\lambda(t-s)} (\lambda + C_1) \|\mathcal{U}(s)u_0\| ds \\ &= \Theta(t) + \int_{t_0}^t e^{-\lambda(t-s)} (\lambda + C_1) \|\mathcal{U}(s)u_0\| ds \\ &\le 1 + \int_{t_0}^t w(t,s) \|\mathcal{U}(s)u_0\| ds, \end{aligned} $$

where

$$ w(t, s) = e^{-\lambda(t-s)}(\lambda + C_1), \quad \text{for } t \ge s \ge t_0. $$

By Lemma 3.1, we obtain that

$$ \begin{aligned} \|\mathcal{U}(t)u_0\| &\le \exp \left( \int_{t_0}^t \left[ (\lambda + C_1) \left[ 1 - \lambda \int_{t_0}^s e^{-\lambda(s-r)} dr \right] \right] ds \right) \\ &= \exp \left( \int_{t_0}^t \left[ (\lambda + C_1) \left[ 1 - e^{-\lambda s} (e^{\lambda s} - e^{\lambda t_0}) \right] \right] ds \right) \\ &= \exp \left( \int_{t_0}^t \left[ (\lambda + C_1) e^{-\lambda(s-t_0)} \right] ds \right) \\ &= \exp \left( \frac{(\lambda + C_1)}{\lambda} (1 - e^{-\lambda(t-t_0)}) \right) \\ &\le \exp \left( \frac{\lambda + C_1}{\lambda} \right). \end{aligned} $$

Consequently, for each $u_0 \in \Lambda_+$, we have

$$ \limsup_{t \to +\infty} \|\mathcal{U}(t)u_0\| \le \tilde{C}_2, \quad \text{where} \quad \tilde{C}_2 = \exp \left( \frac{\lambda + C_1}{\lambda} \right). $$

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

The following is the main result in this Section, and it is an immediate consequence of Proposition 5.1, Theorem 3.2, and Theorem 5.1.

**Theorem 5.2.** There is a global attractor $\mathcal{A}$ corresponding to the semiflow $(\mathcal{U}(t))_{t \ge 0}$. More precisely, $\mathcal{A}$ is a compact, connected, and invariant subset of $\Lambda_+$ that attracts all bounded subsets of $\Lambda_+$.
# 6. Numerical visualization of the global attractor

Visualizing the global attractor of a partial differential equation through numerical simulations can be complex. The complexity comes from several factors. First, PDEs often describe systems with high-dimensional state spaces, requiring a large number of variables to accurately represent the behavior of the system. As a result, numerical simulations need to discretize the state space, leading to a substantial increase in computational resources required as the dimensionality grows. Second, the evolution of the PDEs over time demands solving a system of differential equations numerically on complex geometries, typically using methods such as finite difference, finite element, or graph based methods. These numerical schemes involve computations that scale with the size of the discretization, leading to increased computational costs for larger systems or complex geometries. Additionally, the time span required to observe the long-term behavior of the attractor might be extensive, requiring prolonged simulations. Consequently, exploring the global attractor computationally often necessitates substantial computational resources, efficient numerical algorithms, and parallel computing techniques to tackle the computational complexity and achieve meaningful visualizations.

In order to visualize the global attractor we simulate the model described above on a real sandy loam soil sample captured using micro tomographic imaging. For a comprehensive understanding of the soil samples and the techniques employed to generate CT images, readers are referred to the study conducted by Juyal et al. (2018) [25].

The pore space is initially represented as a 3D binary image with dimensions of 512 × 512 × 512. Figure 2 presents a random $z$ plan of the segmented 3D image; the pore space voxels in the image are identified and labeled as black. The 3D image is presented using a uniform resolution of 24 µm.

In order to address the computational limitations, we extract a 3D portion of the original image measuring 503³ pixels located in the [[50, 100], [50, 100], [150, 200]] region. The porosity of the extracted portion is 0.13%. To simplify the geometric complexity of the problem, we employ a method outlined in [26] to approximate the pore space using spheres. Subsequently, we utilize the numerical model described in [12] to simulate the long-term behavior of the system within the intrinsic geometry of the resulting pore network model. The numerical method achieves efficient computation time and generates highly accurate simulations of the model described earlier in a complex network of geometric primitives. In this discussion, we first provide an overview of the geometric technique employed to extract the minimal set of maximal spheres that cover the pore space. Then, we delve into the specifics of the numerical framework utilized for simulating the global attractor of the problem.

## 6.1. Pore network extraction

To facilitate reader comprehension, we present a summarized description of the geometric modeling of the pore space as a minimal set of balls that cover the pores of the 3D image.
<figure><img src="image_4.png" /><figcaption>Figure 2: Random z-plane slice from the 3D X-ray CT scan of the sandy loam soil sample. The pore space is shown in black. This image represents the structural domain used as the simulation space.</figcaption></figure>

The methodology employed in this study, as detailed in [26], involves using a minimal set of balls to approximate the complex pore space in the 3D image. This approach offers a more concise representation of the pore space, which is better suited for numerical simulations compared to the original set of voxels [12]. It is also considered more realistic than idealized pore network models, as highlighted in [27]. In our research, we adopt balls as the primary primitives for this purpose. However, alternative primitives, such as ellipsoids, could have been incorporated within the same numerical simulation framework [30].

The main principle of the numerical scheme consists in describing the pore space (the domain of the PDE system) by a relational attributed graph in order to discretize the PDE system. Each node of the valuated graph is attached to a geometric primitive (voxel, ball, ellipsoid...) and each arc to an adjacency relationship between two primitives as in [12, 26, 27]. The simplest and direct way is to represent pore space by the graph of the set of voxels as in reference [12] but it leads to high computing time when implementing dynamics simulation. The ball based representation is a good trade-off between the compacity of the model and its precision
[12, 26]. Its basic principle described through details in [12, 26, 27] consists in determining the minimal set of balls (in the cardinal sense) recovering the skeleton of the shape (locus of the centers of all maximal balls, in the inclusion sense, of the shape). Intuitively it is equivalent to inflate soap balls into the main cavities of the shape [12, 26, 27]. If the skeleton is a compact set (in the mathematical sense), by definition of a compact set, it exists at least one finite set of balls recovering the skeleton. Then we describe the shape by the minimal (in the cardinal sense) set(s) of maximal balls recovering the skeleton. We have shown in reference [29] that it exists at least one minimal set recovering the skeleton formed of maximal balls (in inclusion sense). The set of all maximal balls can be computed thanks to the 3D Delaunay triangulation of the boundary points of the shape (pore space). In a first stage we prune the Delaunay tetrahedra by selecting only the tetrahedra included within the shape. The set of the corresponding Delaunay spheres correspond to the set of the maximal balls of the shape. Afterward heuristic algorithms can provide the desired ball network. The principle of the algorithms is to place iteratively the biggest maximal ball such that it intersects no other balls [12, 26, 27]. Practically this simple scheme gives a minimum set of balls, disjoint or tangent two by two, which includes the skeleton of the shape. Each ball captures a local cavity or channel, and the collection of balls collectively approximates the entire pore domain while preserving important morphological features such as volume and connectivity. Using these primitives, we construct an attributed adjacency graph in which each node corresponds to a sphere (representing a pore) and each edge represents adjacency between neighboring spheres. This graph serves as the simulation space, allowing diffusion and microbial processes to be modeled with significantly reduced computational cost.

By considering the collective set of primitives, we obtain an approximation of the overall pore space that enables a comprehensive analysis of connectivity and relationships between different pore-like structures. Figure 3 illustrates the reconstructed pore network corresponding to the selected 3D section, visualized using a Matlab routine.

## 6.2. Numerical schemes of the global model in the extracted pore network

As in [12], we divide the dynamics to diffusion and transformation processes which suit very well for the graph-based representation of the domain. The pore space is presented using a graph of geometric primitives; balls in this experience (see Figure 3).

So given a set $\mathbf{P} = \{p_1, \dots, p_n\}$ of geometric primitives obtained using the approach described above [Olivier Monga, 2007], we construct a graph $G(N,E)$, where: $N = \{1, \dots, n\}$ is the index set of the graph which corresponds to the primitives in $\mathbf{P}$, and $E = \{(i, j)/p_i \cap p_j \neq \emptyset\}$ is the set of edge indices which encodes the geometrical adjacency between the primitives.
<figure><img src="image_5.png" /><figcaption>Figure 3: The minimal set of maximal sphere network covering the pore space of the selected portion from the original 3D binary image.</figcaption></figure>

Let $t \ge 0$, for each node i of N we attach the following vector

$$X_i(t) := \{x_i, y_i, z_i, v_i, b^i(t), n^i(t), m_1^i(t), m_2^i(t), c^i(t)\},$$

where $(x_i, y_i, z_i)$ are the coordinates of the gravitation center of the geometrical primitive $p_i$ and $v_i$ its volume, $b^i(t), n^i(t), m_1^i(t), m_2^i(t), c^i(t)$ are respectively the total mass of MB, DOM, SOM, FOM, CO$_2$ contained within the primitive $p_i$ at time t.

We firstly present the mathematical model for diffusion in the set of geometrical primitives, then we adjust the global model in order to give the general framework for simulating the global mathematical model in the graph of geometrical primitives.

### 6.2.1. Diffusion processes

Fick's first Law of Diffusion describes the rate of diffusion through a medium in terms of the concentration gradient. It states that the flux (J) of particles across a unit area perpendicular to the direction of diffusion is proportional to the concentration gradient.
Specifically, the concentration gradient at a given point in the contact area between two adjacent regular (i.e., geometrically regular) primitives $p_i$ and $p_j$ is given by

$$ \frac{\frac{m_i}{v_i} - \frac{m_j}{v_j}}{L_{i,j}}, $$

where $\frac{m_i}{v_i}$ and $\frac{m_j}{v_j}$ denote the concentrations in primitives $p_i$ and $p_j$, respectively, and $L_{i,j}$ represents the distance between their centers of mass. Since the primitives are regular, the distance between their gravitational centers is used.

Then, the total flux of mass from the primitive $p_j$ to the primitive $p_i$ during a time period $\delta t$ is the following,

$$ J_{i,j} = -DA_{i,j}\delta t \frac{c_i - c_j}{L_{i,j}}, $$

where $A_{i,j}$ is the area of contact between the two primitives, D is the diffusion coefficient of the compound to be diffused and $c_i = \frac{m_i}{v_i}$.

The total variation of mass by diffusion processes at the primitive $p_i$ during a time period $\delta t$ is the sum of all the fluxes coming from the neighboring primitives, which is expressed as follows;

$$ \delta m_i = \sum_{j \in N/(i,j) \in E} -DA_{i,j}\delta t \frac{c_i - c_j}{L_{i,j}}, $$

which results in,

$$ \begin{aligned} \frac{dm_i}{dt}(t) &= \sum_{j \in N/(i,j) \in E} -DA_{i,j} \frac{c_i(t) - c_j(t)}{L_{i,j}} \\ &= D(c_i(t) \sum_{j \in N/(i,j) \in E} -Q_{i,j} + \sum_{j \in N/(i,j) \in E} Q_{i,j}c_j(t)) \\ &= D\left(\frac{m_i(t)}{v_i} \sum_{j \in N/(i,j) \in E} -Q_{i,j} + \sum_{j \in N/(i,j) \in E} Q_{i,j} \frac{m_j(t)}{v_j}\right), \end{aligned} $$

where $Q_{i,j} = \frac{A_{i,j}}{L_{i,j}}$ is the diffusional conductance between primitives $p_i$ and $p_j$.

Gathering the masses $m_i$ of all the primitives in the same vector $M(t) = (m_i(t))_{i \in N}$, we get the following equation,

$$ \frac{dM}{dt}(t) = D\hat{\Delta}M(t), $$

where

$$ \hat{\Delta}_{i,i} = \frac{1}{v_i} \sum_{j \in N/(i,j) \in E} -Q_{i,j}, $$

and

$$ \hat{\Delta}_{i,j} = \frac{Q_{i,j}}{v_j} \text{ for } i \neq j. $$
In practice, the coefficients $A_{i,j}$ and $L_{i,j}$ cannot be determined with high accuracy. Following the approach in [12], we instead compute effective values of $Q_{i,j}$, referred to as diffusional conductance coefficients, which serve as surrogates for the exact geometrical parameters.

### 6.3. Transformation processes and general model

Similarly to the model described above we consider that FOM and SOM are decomposed rapidly and slowly, respectively. DOM comes from the hydrolysis of SOM and FOM. DOM diffuses through water paths (water-filled spheres) and is consumed by MB for its growth. We hypothesized that MB does not move. Dead microorganisms are recycled into SOM and DOM. MB respires by producing CO$_2$.

The changes of the set of biological features ($b^i(t)$, $n^i(t)$, $m_1^i(t)$, $m_2^i(t)$, $c^i(t)$)$_{i \in N}$, due to transformation of different compounds and diffusion processes of DOM in the water-filled primitives, within a time step $\delta t$ are expressed using the following system of equations (7):

$$ \forall i \in N, \left\{ \begin{aligned} \frac{db^i}{dt}(t) &= \frac{Kn^i(t)}{K_b + n^i(t)}b^i(t) - (\eta + \mu)b^i(t), \\ \frac{dn^i}{dt}(t) &= D_n(\hat{\Delta}N(t))_i + \rho\mu b^i(t) - \frac{Kn^i(t)}{K_b + n^i(t)}b^i(t) + c_1m_1^i(t) + c_2m_2^i(t), \\ \frac{dm_1^i}{dt}(t) &= -c_1m_1^i(t) + (1 - \rho)\mu b^i(t), \\ \frac{dm_2^i}{dt}(t) &= -c_2m_2^i(t), \\ \frac{dc^i}{dt}(t) &= \eta b^i(t). \end{aligned} \right. \qquad (7) $$

In the given context, $N(t) = (n^i(t))_{i \in N}$ denotes the vector of the contained DOM mass within all the geometric primitives at time t, $\eta$ represents the relative respiration rate in units of $day^{-1}$, $\mu$ denotes the relative mortality rate in units of $day^{-1}$, $\rho$ signifies the proportion of MB that returns to DOM while the remaining fraction returns to SOM. Furthermore, $c_2$ and $c_1$ correspond to the relative decomposition rates of FOM and SOM respectively, both measured in units of $day^{-1}$. Additionally, $K$ and $K_b$ represent the maximum relative growth rate of MB and the constant of half-saturation of DOM by MB, respectively, both measured in units of $day^{-1}$ and gC (grams of carbon).

The term $D_n(\hat{\Delta}N(t))_i$ represents the mass variation of DOM caused by the exchange of mass between the primitive $p_i$ and all the connected primitives. Here, $D_n$ denotes the molecular diffusion coefficient of DOM in water, measured in units of $cm^2.d^{-1}$.
### 6.4. Aggregation Techniques for Visualizing the Attractors

In the context of visualizing attractors for high dimensional PDEs, it is often necessary to use aggregation techniques to simplify and clarify the data. Aggregation involves summarizing multiple variables or data points into a single value, which can help to reduce the dimensionality of the system and create more effective visualizations.

To better understand our model, we wanted to look at how the biological parameters change over time. We focused on studying the long-term growth of microorganisms in relation to their need for organic matter and the production of CO₂.

To illustrate our results, we introduce the following application $AGG_1$ that for a distribution of $L^2(\Omega)$ returns the integral sum over omega of the distribution, mathematically defined as,

$$AGG_1: \Lambda_+ \to \mathbb{R}_+^5 \\ (\phi_1, \dots, \phi_5) \mapsto (\int_{\Omega} \phi_1, \dots, \int_{\Omega} \phi_5).$$

This function calculates the overall mass of biological parameters in a given sample by integrating their distribution across the domain $\Omega$.

In the numerical method, each primitive $p_i$ is expressed by the contained masses $\{b^i(t), n^i(t), m_1^i(t), m_2^i(t), c^i(t)\}$ within its volume $v_i$ at time t.

So, the total mass of different biological parameters within the pore space is obtained by the following;

$$AGG_1: \quad (l^2(N))^5 \quad \to \quad \mathbb{R}_+^5 \\ [(b^i)_{i \in N}, (n^i)_{i \in N}, (m_1^i)_{i \in N}, (m_2^i)_{i \in N}, (c^i)_{i \in N}] \quad \mapsto \quad [\bar{B}, \bar{N}, \bar{M}_1, \bar{M}_2, \bar{C}],$$

where

$$[\bar{B}, \bar{N}, \bar{M}_1, \bar{M}_2, \bar{C}] := \left[ \sum_{i \in N} b^i, \sum_{i \in N} n^i, \sum_{i \in N} m_1^i, \sum_{i \in N} m_2^i, \sum_{i \in N} c^i \right].$$

Despite the fact that this approach reduces the attractor from the infinite-dimensional space $(l^2(N))^5$ to a subset in $\mathbb{R}^5$, we further aggregate the data to focus solely on the microorganisms growth affected by the total organic matter and the CO₂ occurring in the soil sample. We also, use projection in order to focus on specific parameters. For this purpose we introduce the following applications:

$$AGG_2: \quad \mathbb{R}_+^5 \quad \to \quad \mathbb{R}_+^3 \\ [\bar{B}, \bar{N}, \bar{M}_1, \bar{M}_2, \bar{C}] \quad \mapsto \quad [\bar{B}, \bar{N} + \bar{M}_1 + \bar{M}_2, \bar{C}]$$

and
$$ \mathcal{P}: \quad \begin{array}{ccc} \mathbb{R}_+^5 & \rightarrow & \mathbb{R}_+^3 \\ [\bar{B}, \bar{N}, \bar{M}_1, \bar{M}_2, \bar{C}] & \mapsto & [\bar{B}, \bar{N}, \bar{C}]. \end{array} $$

In the given context, the variable $\bar{B}$ represents the overall population of microorganisms present in the soil sample. Similarly, the quantity $\bar{N}+\bar{M}_1+\bar{M}_2$ corresponds to the combined mass of organic matter encompassing its various forms, including DOM, DOM, FOM, and SOM. Lastly, $\bar{C}$ represents the amount of carbon that has been released through microbial respiration.

### 6.5. Model parameters and initial conditions modeling

The same biological parameters of Arthrobacter sp. 9R as in [12, 31, 32] were utilized. The parameter values assigned were as follows: $\eta$, representing the relative respiration rate, was set to $0.2 \text{ day}^{-1}$; $\mu$, representing the relative mortality rate, was set to $0.5 \text{ day}^{-1}$; $\rho$, indicating the proportion of MB returning to the DOM, was set to 0.55 (with the remaining portion returning to SOM); $c_2$ and $c_1$, representing the relative decomposition rates of FOM and SOM, were respectively set to $0.3 \text{ day}^{-1}$ and $0.01 \text{ day}^{-1}$; $K$, denoting the maximum relative growth rate of MB, was set to $9.6 \text{ day}^{-1}$; $K_b$, signifying the constant of half saturation of DOM, was set to $0.001 \text{ gC.g}^{-1}$ (grams of carbon per gram).

Different scenarios were conducted on the extracted portion of the original 3D image to explore various initial conditions. DOM and MB were distributed in either a heterogeneous or homogeneous manner, randomly. The total masses of DOM and MB in each scenario are summarized in Figure 4.

In each scenario, the mass of DOM in the network of spheres was determined by selecting a value within the range of $[10^{-7}, 9.10^{-4}] \text{ µg voxel}^{-3}$. The DOM mass was distributed among the water-filled spheres, that represent a volume of $1.73 \text{ mm}^3$, either heterogeneously or homogeneously.

To distribute the microorganisms in a heterogeneous manner, following the realistic bacterial model distribution outlined in [33] and used in previous works [12, 32], a random value between 0.05% and 0.15% of the distributed DOM was chosen to determine the mass of microorganisms. These microorganisms were then distributed in the spheres network as patches, with the number and placement of the patches chosen randomly.

We run simulations of the model from the resulting initial conditions during a period of 918 days using an implicit scheme of the model (7). Figure 5 shows simulations during one year of 4 randomly chosen scenarios.
<figure><img src="image_6.png" /><figcaption>Figure 4: In each scenario, the initial masses of MB and DOM are shown, with MB represented in blue and DOM in orange.</figcaption></figure>

## 6.6. Computational Cost and Numerical Validation

The simulations of the modeled scenarios were performed on a standard workstation equipped with an AMD Ryzen 7 PRO 6850H processor and 32 GB of RAM. Although the simulated time span extended over 918 days, the average runtime for completing a full scenario was only about 8 minutes. This highlights a remarkable reduction in computational demand compared to traditional voxel-based simulations [12]. These references also provide detailed validation of the numerical scheme, including convergence, stability, accuracy analyses, and benchmark comparisons, ensuring the reliability of the simulation results.

The gain in efficiency is mainly attributed to the dimensionality reduction achieved by the sphere-based modeling approach. In this framework, the simulation space was reduced from 21,420 voxels to only 218 spatial units (spheres). Each sphere aggregates the information of multiple voxels while preserving essential morphological and topological features of the pore space. As a result, the number of entities over which diffusion and reaction equations must be solved is drastically lowered [12].

This compact representation not only accelerates computation but also enhances scalability, making it feasible to simulate larger domains or longer temporal dynamics. Consequently, the approach enables the study of more complex soil structures without incurring the prohibitive costs typically associated with high-resolution voxel-based simulations.

# 7. Results and discussion

Figures 4 and 5 illustrate the dynamics of microbial biomass (MB), dissolved organic matter (DOM), soil organic matter (SOM), fresh organic matter (FOM), and CO₂
Simulations for 4 Scenarios - Full and Zoomed Views

<figure><img src="image_7.png" /><figcaption>Figure 5: Simulation of some scenarios during one year: Each row corresponds to one scenario, with two subplots: **Left:** full simulation duration of 365 days, **Right:** zoom-in on the first 20 days to highlight initial dynamics. The line colors correspond to the different components as follows: MB (blue), DOM (orange), SOM (green), FOM (red), and CO₂ (purple).</figcaption></figure>

under different scenarios. In Figure 4, the initial masses of MB and DOM for each scenario are shown, highlighting the variability in starting conditions across scenarios. Figure 5 presents the temporal evolution of these components over one year. The left subplots display the full 365-day simulation, while the right subplots focus on the first 20 days to capture the early dynamics.

Comparing scenarios in Figure 5, it is evident that scenarios with higher initial MB or DOM exhibit faster early consumption of DOM and more rapid accumulation of CO₂, indicating higher microbial activity. Conversely, scenarios with lower initial MB show slower DOM uptake and delayed CO₂ production. SOM depletion is more pronounced in scenarios where microbial biomass is initially high, reflecting stronger
decomposition rates. These comparisons underscore the impact of initial conditions
on the temporal progression of carbon transformations in the soil system.

All distributions tend to have the same patterns and to converge towards a
common region and exhibit oscillations within it. The attractor lies in the plane
corresponding to a dead population of microorganisms, indicating that microbial
growth is limited by carbon availability. This observation is consistent with the
experimental findings in [34, 35], which demonstrated carbon limitation in microbial
growth.

Figure 6 shows the long-term system dynamics following the application of ag-
gregation functions $AGG_1$ and $AGG_2$. The panels illustrate the evolution from the
initial state ($t = 0$ h) to intermediate times ($t = 20$ h, $t = 100$ h) and the final
state at $t = 22,000$ h ($\sim 2.51$ years). The initial positions (red stars) and terminal
positions (red circles) indicate how the system evolves over time, highlighting differ-
ences in trajectories under the two aggregation functions. Notably, the comparison
of intermediate frames demonstrates the rate of convergence toward the attractor
and the distinct transient behaviors induced by $AGG_1$ and $AGG_2$.

Figure 7 presents similar dynamics for $AGG_1$ and $P$. The figure highlights how
the system trajectories differ under these aggregation functions, particularly in the
early dynamics (first 100 hours) and in the final steady-state positions. The patterns
of convergence, spread of trajectories, and final positions provide insight into how
each aggregation function influences the long-term system behavior.

Examining Figures 6 and 7, it is clear that the attractor resides within a well-
defined region, distinct from the surrounding data points. These visualizations con-
firm that the attractor is evident in 3D space.

Figure 8 illustrates the stability of a distribution belonging to the attractor
by showing the final state of a random scenario. Our analysis demonstrates that
when simulations start from this distribution, they consistently remain within the
attractor. This highlights the attractor’s role in governing the system’s long-term
behavior, ensuring that trajectories originating from within it stay tightly clustered
around its basin of attraction.

Despite variations among different scenarios, the overall trend shows convergence
towards a shared region, with periodic fluctuations within this range. This region
aligns with the $\bar{B} = 0$ plane, indicating that microbial growth is restricted due to a
lack of accessible organic matter resources.

8. Conclusion

In conclusion, this paper analyzed a reaction-diffusion model describing micro-
bial decomposition of organic matter in a three-dimensional soil structure. The
model was formulated using nonlinear parabolic PDEs, for which we established
existence and uniqueness of solutions. Furthermore, we demonstrated the presence
of a global attractor.

To simulate the model and visualize its long-term behavior, we provided a nu-
merical tool based on modeling the pore space using an optimal sphere network and
<figure><img src="image_8.png" /><figcaption>Figure 6: Long-term system dynamics following the application of aggregation functions $AGG_1$ and $AGG_2$. Panel (A) shows the initial state at $t = 0$ h (initial total masses). Panel (B) corresponds to $t = 20$ h, (C) to $t = 100$ h, and (D) to $t = 22,000$ h ≈ 2.51 years. Initial positions are marked with red stars and terminal positions with red circles. A concise legend is displayed next to the pictures to improve readability. An animation of all simulation frames is available at this link.</figcaption></figure>

formulated a numerical scheme that corresponds to the model in the sphere network. This tool is suitable for simulating the model over long periods as it requires reasonable computing power and time.

We provided aggregation functions that reduce the dimensionality of the solution to better visualize the global attractor of the microbial system in soil. The simulations presented in this paper were conducted on a fully saturated real sandy loam soil sample captured using a micro-CT scanner and segmented in the form of a binary 3D image.

We showed that the solution of the model, regardless of the initial distribution, converges to the attractor, which is a region where microorganisms die due to the absence of organic matter.

The main contribution of the present work is to present a mathematical analysis of the PDE system ruling the microbial decomposition in porous media. The principal outcome is that this system admits a global attractor. This PDE system includes diffusion processes (Laplacian) and transformation processes (biological transformations) in a complex space (pore space). The complex shape describing the pore space is represented by a set of voxels determined from Computed Tomography images. In a mathematical point of view, these voxels can be considered as cubes in 3D space. Practically the number of voxels describing pore space is up to dozens and hundred
<figure><img src="image_9.png" /><figcaption>Figure 7: Long-term system dynamics following the application of aggregation functions $AGG_1$ and $P$. Panel (A) shows the initial state at $t = 0$ h (initial total masses). Panel (B) corresponds to $t = 20$ h, (C) to $t = 100$ h, and (D) to $t = 22,000$ h $\approx 2.51$ years ($\sim 2$ years, 6 months, 4 days). Initial positions are marked with red stars; red circles indicate terminal positions in each frame. An animation of all frames is available at this link.</figcaption></figure>

of millions that makes difficult and even impossible the direct solving of the PDE system using existing solvers (Freem++,...). As far as we know, very few works deal with a rigorous mathematical analysis of this PDE system. In the reference [5], one gives argument for proving the convergence of the system. In the present work, we showed, using the semi-group theory, that the system admits a global attractor. We used a discretized version of the PDE system by means of forward and backward Euler schemes running on a valuated graph representing the dynamics space in a compact way. We present simulation results on hundreds of data sets to illustrate the convergence towards an attractor.

In terms of applications, our work can inform soil management and carbon cycling studies, helping to predict microbial responses under varying conditions. For future perspectives, this framework could be extended to multi-phase flow systems, different soil textures, or coupled nutrient-microbial interactions. Furthermore, it can support the development of predictive models for soil carbon sequestration and optimization of soil amendment strategies.

**Perspectives:** One motivating perspective is to find out cases where the system is vibrating. Practically that would mean that the CO2 coming from the breathing would converge to an asymptote in order to avoid the complete removing of the
<figure><img src="image_10.png" /><figcaption>Figure 8: Long-term system dynamics following the application of the aggregation functions $AGG_1$ and $P$. All initial states are marked with stars, and their corresponding terminal states with circles. Red stars denote randomly chosen initial distributions, while the black star represents an initial distribution taken from the final states of simulations in Figures 6 and 7 (i.e., a distribution lying in the proposed global attractor). The black circle indicates its terminal state, which coincides with the initial point, showing that this distribution remains in the attractor.</figcaption></figure>

biomass. In the case where no vibrating state could be reached, that would mean a simple convergence of the system for each case. Another perspective would be to analyze the solution of the system for any time and not only for long time (equivalent to infinite).

## 9. Appendix: Discretization, Numerical Implementation of the PDE System, and Time Convergence

In order to discretize the PDE system (1), we used forward and backward Euler discretization scheme. The discretization of the diffusion process is realized thanks to first Fick law. In our practical implementation same as in reference [25], we assumed that only the dissolved organic matter (DOM) diffuse and therefore that $D_n = D_c = 0$. The PDE system becomes:

$$ 
\left\{ 
\begin{aligned}
\frac{\partial b}{\partial t} &= -\eta b - \mu b + \frac{K_n}{K_b+n} b \\
\frac{\partial n}{\partial t} &= D_c \Delta(n) + \rho \mu b - \frac{K_n}{K_b+n} b + c_1 m_1 + c_2 m_2 \\
\frac{dm_1}{dt} &= (1-\rho)\mu b - c_1 m_1 \\
\frac{dm_2}{dt} &= -c_2 m_2 \\
\frac{dc}{dt} &= \eta b
\end{aligned}
\right.
\quad (A-2) $$

In order to discretize the equation (A-2) we split the system into a transformation part (system (A-3)) and a diffusion part (system (A-4)). We got the following two systems:
$$ \left\{ \begin{array}{l} \frac{\partial b}{\partial t} = -\eta b - \mu b + \frac{Kn}{K_b+n} b \\ \frac{\partial n}{\partial t} = \rho \mu b - \frac{Kn}{K_b+n} b + c_1 m_1 + c_2 m_2 \\ \frac{dm_1}{dt} = (1-\rho)\mu b - c_1 m_1 \\ \frac{dm_2}{dt} = -c_2 m_2 \\ \frac{dc}{dt} = \eta b \end{array} \right. \qquad (A-3) $$

$$ \frac{\partial n}{\partial t} = D_c \Delta(n) \qquad (A-4) $$

<figure><img src="image_11.png" /><figcaption>Figure 9: Perspective views of the ball network (1526 balls).</figcaption></figure>

We discretize in time using Euler discretization schemes the PDE systems (A-3) and (A-4). For system (A-3) we used Euler forward discretization scheme and for system (A-4) Euler backward discretization scheme. The reason why using Euler backward discretization scheme for the diffusion is that forward discretization scheme would imply too little time steps in order to avoid negative values involving very high computing times (see reference [25]).

In a first step we will recall the principles of forward and backward Euler discretization scheme (see MIT courses: link).

In a second step we applied Euler forward and backward schemes respectively to subsystems (A-2) and (A-3).

Given a PDE: $\frac{\partial(u(x,y,z,t))}{\partial t} = F(u, x, y, z, t)$, the forward Euler scheme is defined by: $u(t + \Delta t, x, y, z) = u(t, x, y, z) + \Delta t F(u(t, x, y, z), x, y, z, t)$. Therefore the principle is to compute directly the solution at time $t + \Delta t$ from the solution at time t using a first order approximation (first order Taylor series expansion). The approximation
error is proportional to the time step $\Delta t$. The stability limitation is that it may require very small time steps for stability especially for stiff problems.

The backward Euler method is defined by: $u(t + \Delta t, x, y, z) = u(t, x, y, z) + \Delta t F(u(t + \Delta t, x, y, z), x, y, z, t + \Delta t)$. Therefore the principle is to compute the solution at time $t$ from solution at time $t + \Delta t$. It requires solving a system of equations to find the solution at next time step. Same as for forward Euler scheme the approximation error is proportional to time step $\Delta t$. For stiff problems, it allows to use bigger time steps than the forward scheme. For linear problems the backward scheme is often unconditionally stable.

The discretization of PDE system (A-2) using Euler forward scheme is:

$$ 
\left\{ 
\begin{array}{l}
b(t + \delta t) = b(t) - \eta b(t) \delta t - \mu b(t) \delta t + \frac{Kn(t)}{\kappa_b + n(t)} b(t) \delta t \\
n(t + \delta t) = n(t) + \rho \mu b(t) \delta t - \frac{Kn(t)}{\kappa_b + n(t)} b(t) \delta t + c_1 m_1(t) \delta t + c_2 m_2(t) \delta t \\
m_1(t + \delta t) = m_1(t) + (1 - \rho) \mu b(t) \delta t - c_1 m_1(t) \delta t \\
m_2(t + \delta t) = m_2(t) - c_2 m_2(t) \delta t \\
c(t + \delta t) = c(t) + \eta b(t) \delta t
\end{array}
\right. 
\quad (A-5) 
$$

The discretization of equation (A-4) is derived from first Fick law defining the flow between two points along a diffusion process. The first Fick's law states that the diffusion flow of a substance from two adjacent points is directly proportional to the difference of the substance concentrations between the two points. That means that the substance move from areas of higher concentrations to areas of lower concentrations with a speed proportional to the concentration gradients. First Fick law can be expressed as (see reference [28]):

$$ J = -D_c \frac{\partial n}{\partial x}, $$

where, $J$ represents the diffusion flux (amount of substance per unit area per unit time), $D_c$ is the diffusion coefficient of the substance, and $\frac{\partial n}{\partial x}$ is the concentration gradient along the diffusion direction. The use of first Fick law for defining an implicit Euler numerical scheme for implementing diffusion processes is described precisely in reference [25]. In the present paper we recall the practical numerical scheme we get and also show numerical simulations showing the convergence in time of the whole process.

As described precisely in references [25, 12], we described the pore space by means of an adjacency relational graph (balls network) where each node is attached to a ball. Afterward, the solving of PDE system (A-2) is reduced to graph updating. We used a splitting scheme for transport (diffusion, system (A-4)) and reaction (transformation, system (A-3)).

We expressed Euler backward numerical scheme by this way. We set: $\Theta_{ij} = D_c \frac{s_{ij}}{d_{ij}} \delta t$ where $D_c, s_{ij}, d_{ij}, \delta t$ are respectively the diffusion coefficient, the area of the contact surface between balls i and j, the distance between the two balls centers i and j, and $\delta t$ the discretization time step.
<figure><img src="image_12.png" /><figcaption>Figure 10: Kinetics of microbial decomposition for various discretization time steps $\delta t$ (sequential transformation and diffusion) over a 5-day simulation. Panels are ordered from left to right and top to bottom for decreasing $\delta t$: 3600 s, 1800 s, 900 s, 450 s, 225 s, 112.5 s, 56.25 s, 28.125 s, 14.0625 s, 7.03 s, 3.51 s, and 1.76 s. The X-axis shows time (hours) and the Y-axis shows the percentage of initial mass. Curves: DOM (green), MB (blue), CO$_2$ (red), SOM (light blue). Normalized Cross Correlation (NCC) values between successive CO$_2$(t) curves indicate convergence, reached when NCC is within $10^{-6}$ of 1 for four consecutive time steps.</figcaption></figure>

In order to define automatically the time step, we used the following principle. We started with a high time step and divided it iteratively by two until the convergence was reached. The convergence is tested by calculating the normalized intercorrelation between two successive CO$_2$(t) curves.

Figure 10 presents the results obtained using the same dataset as the one of Figures 10–13. Figure 13 presents the same results as Figure 10 but where each kinetic (MB, DOM, CO$_2$, SOM) is displayed for the different time steps within the same frame axis.
We started with a time step set to 3600 seconds and decreased it iteratively by dividing it by two (3600 s, 1800 s, 900 s, 450 s, ...). For each time step, we calculated the Normalized Cross Correlation (NCC) of the two CO₂(t) curves corresponding to the time step and to the previous one. We assumed that convergence was reached when the NCC reached 1 ± 10⁻⁶ for four successive time steps.

We found that the convergence is reached for a time step of about 30 s, which also confirms the results reported in [10].

<figure><img src="image_13.png" /><figcaption>Figure 11: Temporal convergence analysis. L² error with respect to the finest time-step ($\delta t$ = 0.439s) solution for (MB), (DOM), (SOM), (CO₂), and their total sum, as a function of the physical time step $\delta t$.</figcaption></figure>

To explicitly assess the temporal convergence of the numerical scheme, the solution obtained with the finest temporal resolution is used as a reference. For each coarser time step, the numerical error is quantified by computing the L² norm of the difference between the solution and the reference solution over the entire simulation time. This error is evaluated separately for the four state variables: (MB), (DOM), (SOM), and (CO₂), as well as for their total sum.

Figure 11 presents the resulting temporal convergence curves. Error values are plotted at equally spaced positions for clarity, while the x-axis labels indicate the
exact physical time step $\delta t$ (in seconds). The results show a consistent decrease of the error as the time step is refined, providing clear quantitative evidence of temporal convergence of the proposed numerical approach.

<figure><img src="image_14.png" /><figcaption>(a) Simulation over 30 days with $\delta t = 32.4$ s and a relative SOM decomposition rate set to 0.</figcaption></figure>

<figure><img src="image_15.png" /><figcaption>(b) Simulation over 300 days with $\delta t = 32.4$ s and a relative SOM decomposition rate set to 0.01.</figcaption></figure>

Figure 12: Long-term simulations using a sufficiently small time step for temporal convergence ($\delta t = 32.4$ s). X-axis: time in hours; Y-axis: percentage of total initial mass. In both cases, the solution converges toward a stable asymptotic state. The two panels illustrate the impact of extending the simulation time horizon and modifying the relative rate of SOM decomposition. The computational cost increases significantly for longer simulations, with CPU times of 2.25 minutes for the 30-day simulation and 70 minutes for the 300-day simulation on a regular PC (MacBook Pro).

An important practical consequence of this temporal convergence analysis is that, once convergence has been established, larger physical time steps can be safely employed to extend the simulations over much longer periods without significant loss of accuracy. As shown in Figure 11, the numerical error rapidly decreases as the time step $\delta t$ is refined, and remains bounded for time steps larger than the convergence threshold. This property allows a compromise between accuracy and computational efficiency. As an illustration, Figure 12 illustrates the capability of the numerical scheme to perform long-term simulations once temporal convergence has been established. The left panel shows a simulation over 30 days, while the right panel extends the simulation to 300 days using the same convergent time step $\delta t = 32.4$ s. Despite the substantially longer time horizon and increased computational cost, the numerical solution remains stable and converges toward a physically consistent asymptotic state. These results demonstrate that, after convergence is achieved, the use of an appropriately chosen time step allows prolonged simulations without compromising numerical accuracy, which is essential for studying long-term biogeochemical dynamics.

If we denote by $c_i$, $v_i$, $N_i$, $\vartheta(N_i)$, respectively, the concentration at node $i$, the volume of the node $i$, the node $i$, and $\vartheta(N_i)$, the neighbor of the node $N_i$, the relationship between concentrations at successive iterations (time step $\delta t$) can be expressed as follows:
$$ \begin{bmatrix} c_1 \\ c_2 \\ \vdots \\ c_p \end{bmatrix}^k = \begin{pmatrix} \frac{1}{v_1} & \cdots & 0 \\ \vdots & \ddots & \vdots \\ 0 & \cdots & \frac{1}{v_p} \end{pmatrix} \begin{pmatrix} v_1 + \sum_{N_j \in \vartheta(N_1)} \Theta_{1,j} & \cdots & -\Theta_{1,p} \\ \vdots & \ddots & \vdots \\ -\Theta_{p,1} & \cdots & v_p + \sum_{N_j \in \vartheta(N_p)} \Theta_{p,j} \end{pmatrix} \begin{bmatrix} c_1 \\ c_2 \\ \vdots \\ c_p \end{bmatrix}^{k+1} \quad (A-6) $$

where $c_i^k = (1 + \frac{1}{v_i} \sum_{N_j \in \vartheta(N_i)} \Theta_{ij}) c_i^{k+1} - \frac{1}{v_i} \sum_{N_j \in \vartheta(N_i)} \Theta_{ij} c_j^{k+1}$, p is the number of regions (graph nodes), k is the iteration number ($t = t_0 + k\delta t$; $t_0$ is the initial time).

<figure><img src="image_16.png" /><figcaption>Figure 13: Same as figure 10 but where MB, DOM, CO2, SOM are represented in the same frame for each discretization time step.</figcaption></figure>

We solved the above system using Conjugate Gradient Method (CGM). We stress that the equation system (A-6) refers to the diffusive step only because the microbial dynamic step cannot be expressed directly thanks to implicit numerical scheme.

In the present paper versus previous ones [31, 10, 12], we ran complementary simulation showing the convergence in time of the whole process (diffusion, transformation) for the ball based pore space geometric model.

In order to illustrate the convergence in time of the whole process we extracted from the dataset presented in reference [25] a sub image of sizes 150x150x150 of the
whole image of sizes 512x512x512. The pore space of the sub image is represented by a ball network including 1526 balls (see figure 9).

<figure><img src="image_17.png" /><figcaption>Figure 14: Kinetics for different spatial resolutions of the ball network (24 µm, 48 µm, 96 µm, 192 µm). Same format as Figure 14. Panels are ordered from left to right and top to bottom with spatial resolutions 24 µm, 48 µm, 96 µm, 192 µm. Same format as Figure 14 but for different spatial resolutions of the ball network. From left to right and top to bottom, the spatial resolutions are respectively: 24 µm, 48 µm, 96 µm, 192 µm.</figcaption></figure>

In order to illustrate the spatial convergence, we ran simulations using the ball network at different spatial resolutions. We took the same initial distribution of DOM and MB: DOM is spread homogeneously in the ball network and MB is placed in the largest ball. We considered the ball network at different resolutions by iteratively multiplying the resolution by two (equivalent to dividing the balls' coordinates and radius by two). For each resolution, we updated the diffusion coefficient and the half-saturation constant $K_b$, expressed using voxels as the space unit. We ran simulations using the ball network at scales: 24 µm (initial resolution), 48 µm, 96 µm, and 192 µm. The results clearly illustrate the spatial convergence of the simulation process and also its reasonable stability with respect to the spatial resolution. Figures 14 and 15 show the simulation results using the same display principle as in previous Figures 10 and 13.

We put 15µg of Dissolved Organic Matter (DOM) and 0.045µg of biomass (MB) into the ball network. DOM was put homogeneously and MB into the biggest ball
<figure><img src="image_18.png" /><figcaption>Figure 15: Dynamics comparisons for different spatial resolutions (24 µm, 48 µm, 96 µm, 192 µm) for the different kinetics (MB, DOM, SOM, CO₂). Same format as Figure 13 but for different spatial resolutions.</figcaption></figure>

(radius 10.37). In order to simulate the microbial mineralization of organic matter, we took the biological parameters considered in [10, 31, 12] for *Arthrobacter* sp. 9R. The relative respiration rate η was set to 0.2 day⁻¹. The relative mortality rate μ was set 0.5 day⁻¹. The proportion of MB returning to DOM ρ was set to 0.55 (with the remaining portion returning to SOM). The maximum relative growth rate of MB K was set to 9.6 day⁻¹. The constant of half-saturation of DOM Kₕ was set to 0.001 gC g⁻¹. c₁ and c₂ the relative decomposition rates of SOM and FOM were respectively set to 0.01 day⁻¹ and 0.3 day⁻¹.

We ran the whole simulation process for decreasing time steps δt in order to illustrate the convergence in time of the scheme. Figure 10 shows the curves describing the dynamics we got for five days period. Figure 12 shows the curves for longer periods when the time step was set to a convergence value.

CRediT authors contribution statement

Mohammed ELGHANDOURI: Conceptualization, mathematical study, mathematical proofs, manuscript writing, review and editing.

Mouad KLAI: Conceptualization, Numerical methodology, formulation and imple-
mentation, scenarios modeling, simulations, manuscript writing, review and editing.

Khalil EZZINBI: Conceptualization, mathematical methodology and supervision.

Olivier MONGA: Conceptualization, numerical methodology, supervision, and manuscript review and editing.

## Acknowledgments

The research described in this article was made possible thanks to the French institute for development (IRD) through the PhD grants offered to Elghandouri and Klai and thanks to the Moroccan CNRST through the project I-Maroc.

We are grateful to Dr. Frédéric Hecht for engaging in insightful discussions regarding the attractor visualization and the temporal convergence analysis.

A special thank you to the Editor-in-Chief and the anonymous reviewers for their insightful comments and valuable suggestions that have greatly improved this work.

## References

[1] L. Philippot, J. M. Raaijmakers, P. Lemanceau, and W. H. Van Der Putten, "Going back to the roots: the microbial ecology of the rhizosphere," Nature Reviews Microbiology, vol. 11, no. 11, pp. 789-799, Nov. 2013.

[2] R. D. Bardgett and W. H. Van Der Putten, "Belowground biodiversity and ecosystem functioning," Nature, vol. 515, no. 7528, pp. 505-511, Nov. 2014.

[3] N. Fierer, "Embracing the unknown: disentangling the complexities of the soil microbiome," Nature Reviews Microbiology, vol. 15, no. 10, pp. 579-590, Oct. 2017.

[4] P. C. Baveye, W. Otten, A. Kravchenko, M. Balseiro-Romero, È. Beckers, M. Chalhoub, C. Darnault, T. Eickhorst, P. Garnier, S. Hapca, and S. Kiranyaz, "Emergent properties of microbial activity in heterogeneous soil microenvironments: different research approaches are slowly converging, yet major challenges remain," Frontiers in Microbiology, vol. 9, Art. no. 367238, Aug. 27, 2018.

[5] B. Lèye, C. Goudjo, et M. Sy. Convergence analysis of a parabolic nonlinear system arising in biology. Afrika Matematika, 2013, vol. 24, no 2, p. 179-194.

[6] C. M. Dafermos. Semiflows associated with compact and uniform processes. Mathematical systems theory, 1974, vol. 8, no 2, p. 142-149.

[7] Genty A, Pot V. Numerical simulation of 3D liquid-gas distribution in porous media by a two-phase TRT lattice Boltzmann method. Transport in porous media. 2013 Jan;96:271-94.

[8] Haddad, W. M. et Chellaboina, V. Stability and dissipativity theory for non-negative dynamical systems: a unified analysis framework for biological and physiological systems. Nonlinear Analysis: Real World Applications, 2005, vol. 6, no 1, p. 35-65.
[9] Monga O, Bousso M, Garnier P, Pot V. 3D geometric structures and biological activity: Application to microbial soil organic matter decomposition in pore space. *Ecological Modelling*. 2008 Sep 10;216(3-4):291-302.

[10] Monga O, Hecht F, Serge M, Klai M, Bruno M, Dias J, Garnier P, Pot V. Generic tool for numerical simulation of transformation-diffusion processes in complex volume geometric shapes: application to microbial decomposition of organic matter. *Computers & Geosciences*. 2022 Dec 1;169:105240.

[11] Belghali Z, Monga O, Klai M, Abdelwahed EH, Druoton L, Pot V, Baveye PC. Geometric modelling of 3D pore space using curve skeleton: Application to computational microbiology of soil organic matter mineralization. *PLoS One*. 2025 Nov 7;20(11):e0331031.

[12] Klai M, Monga O, Jouini MS, Pot V. A voxel-based approach for simulating microbial decomposition in soil: Comparison with LBM and improvement of morphological models. *PloS one*. 2025 Mar 3;20(3):e0313853.

[13] Klai M. Simulation and Analysis of Joint Transformation and Diffusion Dynamics in Irregular Domains: Application to Computational Modeling of Microbial Decomposition in Porous Media (Doctoral dissertation, Sorbonne Université; Université Cadi Ayyad (Marrakech, Maroc)).

[14] Lèye B, Nguyen-Ngoc D, Monga O, Garnier P, Nunan N. Simulating biological dynamics using partial differential equations: application to decomposition of organic matter in 3D soil structure. *Vietnam Journal of Mathematics*. 2015 Dec;43:801-17.

[15] D. D. Bainov and P. S. Simeonov. Integral inequalities and applications, vol. 57, Springer Science & Business Media, 2013.

[16] A. Batkai, M. K. Fijavz, et A. Rhandi. Positive operator semigroups. Operator Theory: advances and applications, 2017, vol. 257.

[17] Hale JK. Asymptotic behavior of dissipative systems. American Mathematical Soc.; 2010 Jan 4. (1988).

[18] J. Milnor. On the concept of attractor. The theory of chaotic attractors. Springer, New York, NY, 1985. 243-264.

[19] D. Nguyen-Ngoc et al. Modeling microbial decomposition in real 3d soil structures using partial differential equations. *International Journal of geosciences* 2013 (2013).

[20] J. Ning. The global attractor for the solutions to the 3D viscous primitive equations. *Discrete & Continuous Dynamical Systems* 17.1 (2007): 159.
[21] A. Pazy. Semigroups of linear operators and applications to partial differential equations, vol. 44, Springer Science & Business Media, 2012.

[22] V. A. Pliss. Nonlocal problems of the theory of oscillations. Academic Press, 1966.

[23] G. R. Sell, and Y. You. Dynamics of evolutionary equations. Vol. 143. New York: Springer, 2002.

[24] T. Yoshizawa. Stability Theory by Liapunov's Second Method. Vol. 9. Mathematical Society of Japan, 1966.

[25] A. Juyal, T. Eickhorst, R. Falconer, P. C. Baveye, A. Spiers, W. Otten. Control of pore geometry in soil microcosms and its effect on the growth and spread of Pseudomonas and Bacillus sp., Frontiers in Environmental Science, Volume 6, 2018, Article 73.

[26] O. Monga, F. N. Ngom, J. F. Delerue. Representing geometric structures in 3D tomography soil images: Application to pore-space modeling, Computers & Geosciences, Volume 33, Issue 9, 2007.

[27] T. Bultreys, L. V. Hoorebeke, V. Cnudde. Multi-scale, micro-computed tomography-based pore network models to simulate drainage in heterogeneous rocks, Advances in Water Resources, Volume 78, 2015.

[28] A. V. Fick. On liquid diffusion. The London, Edinburgh, and Dublin Philosophical Magazine and Journal of Science, 1855, vol. 10, no 63, p. 30-39.

[29] N. F. Ngom, O. Monga et al. 3D shape extraction segmentation and representation of soil microstructures using generalized cylinders. Computers & geosciences, 2012, vol. 39, p. 50-63

[30] A. T. Kemgue, O. Monga, S. Moto, V. Pot, P. Garnier, P. C. Baveye, A. Bouras. From spheres to ellipsoids: Speeding up considerably the morphological modeling of pore space and water retention in soils, Computers & Geosciences, Volume 123, 2019.

[31] O. Monga, P. Garnier, V. Pot, E. Coucheney, N. Nunan, W. Otten, and C. Chenu. Simulating microbial degradation of organic matter in a simple porous system using the 3-D diffusion-based model MOSAIC. Biogeosciences, 2014, vol. 11, no 8, p. 2201-2209.

[32] B. Mbe, O. Monga, V. Pot, W. Otten, F. Hecht, X. Raynaud, N. Nunan, C. Chenu, P. C. Baveye, P. Garnier. Scenario modelling of carbon mineralization in 3D soil architecture at the microscale: Toward an accessibility coefficient of organic matter for bacteria. European Journal of Soil Science, 2022, vol. 73, no 1, p. e13144.
[33] X. Raynaud, N. Nunan. Spatial ecology of bacteria at the microscale in soil. PloS one, 2014, vol. 9, no 1, p. e87217.

[34] J. L. Soong, L. Fuchslueger, S. Maranon-Jimenez, et al. Microbial carbon limitation: The need for integrating microorganisms into our understanding of ecosystem carbon cycling. Glob Change Biol. 2020; 26: 1953-1961.

[35] T. W. Walker, C. Kaiser, F. Strasser, et al. Microbial temperature sensitivity and biomass change explain soil carbon loss with warming. Nature climate change, 2018, vol. 8, no 10, p. 885-889.