A voxel-based approach for simulating microbial decomposition in soil: Comparison with LBM and improvement of morphological models
Microbial activity in soil is essential for maintaining soil health and supporting ecosystem functioning. Microorganisms contribute to nutrient cycling, organic matter decomposition, and plant growth promotion . However, soil is a highly diverse and dynamic environment, making it challenging to capture its full complexity. Additionally, as a heterogeneous and intricate medium, it poses difficulties for directly observing and measuring microbial activity . Traditional approaches, such as culturing and microscopy, are time-consuming and provide limited insights. Therefore, innovative and advanced numerical simulation techniques are needed to study soil microbial activity . Computed Tomography (CT) imagery offers a non-invasive method for visualizing porous media structures at the microscale of the soil microbial habitats. Standard 3D CT imaging generates large image data volumes, often around one billion voxels . Over the past two decades, significant research has focused on improving computational models for simulating dynamics from 3D CT images of soil. These works aimed to better understand the movement and interaction of soil components in porous and fractured media. These studies primarily integrate transport mechanisms with advanced reaction (transformation) processes . Several methodologies have been proposed to simulate chemical transport phenomena in soil, including the Lattice-Boltzmann Method (LBM) , smooth particle hydrodynamics , hybrid Lattice Boltzmann-direct numerical simulation (DNS) , and pore network geometrical models (PNGM) . Among these, LBM-based approaches have gained popularity for simulating diffusion in porous media using micro 3D CT soil images . However, LBM-based models are computationally expensive due to their particle-based process of collision and propagation, leading to high memory requirements and huge computing times. For example, on a regular laptop, a five-day real-time simulation of microbial decomposition of organic matter with an image size of 512 × 512 × 512 voxels and a porosity of 17% takes approximately three weeks . In contrast, PNGM-based simulations offer a cost-effective alternative with reduced computing time and moderate memory requirements. These models represent pore space as a graph of connected primitives, such as spheres, ellipsoids, or cylinders . However, PNGM simulations require calibration, particularly in optimizing transport (diffusion) by adjusting diffusional conductance coefficients between connected primitives . Despite their computational efficiency, PNGM models are limited by their dependence on such calibrations and may lack precision due to the difference in size distribution of the geometrical primitives used. Zech et al. (2022) simulated microbial decomposition in porous media by coupling diffusion with a microbial decomposition model within a partial differential equation framework, solved using the Galerkin method. Similar to our approach, they numerically solved a PDE-based reaction-diffusion system within pore spaces derived from CT imaging . However, this approach is limited to 2D due to the high computational cost of extending it to 3D. This highlights the ongoing need for scalable 3D solutions that combine the accuracy of voxel-based models with efficient computational strategies. This study addresses key limitations in existing methods by introducing a new computational approach for simulating microbial activity in soil using 3D CT images of pore space. Unlike traditional graph-based techniques modeling pore space with geometrical primitives, our method directly represents the pore space as a graph of connected voxels and models transport phenomena using Fick’s law of diffusion . The resulting improved balance obtained between accuracy and complexity makes it possible to generate data for learning diffusional conductance coefficients in pore network geometrical models. In the Materials and Methods section, we describe the coupling of diffusion processes with a validated model for microbial activity, based on the SOMKO  and MIOR  models, where microbes undergo survival, maintenance, and growth. We detail the used Euler discretization schemes (Sections 2.1-2.4). We then present the calculation of diffusional conductance coefficients in pore network geometrical models (Section 2.5). In Section 2.6 we describe the data set used to compare our new VGA to LBioS and PNGM. We present and discuss the results in Section 3. 

2.1 Simulation principles In this study, we only considered the condition that soil is fully saturated, i.e., soil pore space is filled with water. We model the dynamics of microbial decomposition using five key compounds, as illustrated in : Microbial biomass (MB): Represents the mass of microorganisms in the sample. Microbial respiration (CO2): Indicates the produced carbon dioxide through microbial decomposition, indicative of microbial growth. 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 the hydrolysis of FOM and SOM and biomass recycling, available for microbial uptake or transport within the pore space. Note again that for environments where saturation is not complete, the same simulation framework can be used after a method of drainage  which refers to the method used to redistribute liquid and air within a porous medium, typically by gradually removing liquid from the system. This process mimics the movement of liquid due to capillary forces and gravity, leading to the formation of air-filled pores. In simulations, this is often achieved by progressively reducing the liquid content and adjusting the pore space accordingly. This process allows the model to represent partially saturated conditions before transforming the system into a two-phase configuration by treating the air phase as a non-interacting solid phase . Let I i , j , k i , j , k be a 3D binary image where the voxels forming pore space (void voxels) are tagged by 0 and the voxels attached to solid matter (solid voxels) are tagged by 1. Let V = i , j , k : I i , j , k = 0 = v 1 , … , v n be the set of pore space voxels and N = 1 , … , n the index set of V , where n = c a r d V represents the number of voxels of the pore space. We construct an adjacency valuated graph G V , E where V is the set of nodes, and E = { i , j ∈ N 2 | v i ∩ v j ≠ ∅ } is the set of edges. In this work, we use 6-connecvity between voxels . Let t be the variable representing time. To each voxel v i of the pore space, we attach: the mass of MB ( x i 1 t , the mass of DOM ( x i 2 t , the mass of SOM x || i 3 t , the mass of FOM ( x i 4 t , the mass of CO 2 x || i 5 t . To each node v i let us denote by X i t = x i 1 t , x i 2 t , x i 3 t , x i 4 t , x i 5 t the masses of the compounds contained within voxel v i at time t. We simulate microbial activity by updating the constructed graph according to transformation and diffusion processes described in . The principle is to break down the complex processes, that can be mathematically modelled, into simpler steps. Let us assume that we get a representation of the pore space and its compounds at time t in the form of G t V , E , X i t i ∈ N . To get G t + δ t V , E , X i t + δ t i ∈ N at time t + δ t , after microbial activity processes we apply a transformation model that encodes the biological conversion of different compounds within each voxel. Next, we employ a discretized model to mimic the diffusion of compounds between the connected voxels. The modeling approach is based on dynamics splitting approach, where the diffusion and reaction processes are handled separately. The spatial discretization is fixed according to the computed tomography image resolution. The same framework can be used to model diffusion-transformation processes at any resolution, as long as the parameters are converted based on the voxel size, which depends on the resolution of the 3D image. 2.2 Transport modelling and graph diffusion equation 2.2.1 Fick’s law of diffusion. Fick’s law of diffusion describes the process of molecular or particle movement through a medium, such as gases, liquids, or solids. It is expressed as follows : J = − D d C d x (1) Where J represents the diffusion flux (amount of substance per unit area per unit time), D is the of the substance and d c / d x the concentration gradient along the diffusion direction. 2.2.2 Diffusion modelling: graph diffusion equation. We consider the diffusion of dissolved organic matter (DOM). Let m i t = x i 2 t represent the mass of DOM at time t in the voxel p i . The diffusion process in our framework is modeled using a graph-based approach where voxels represent graph nodes, and edges connect neighboring voxels. The concentration at a voxel corresponds to its mass since the voxel volume is normalized. The mass flow between time t and time t + δ t from the voxel at position p i to its neighboring voxel at position p i is given by the first Fick’s law as follows: J i , j = − D D O M . m i t − m j t . δ t (2) where D D O M is the DOM diffusion coefficient in water. Thus, the DOM mass at voxel i at time t + δ t due to mass exchange with neighbouring voxels is expressed as follows: m i t + δ t − m i t = ∑ j : i , j ∈ E − D D O M δ m i , j t δ t (3) Then the governing equation is: d m i t d t = ∑ j : i , j ∈ E − D D O M m i t − m j t (4) Gathering all the masses of all the voxels of the pore space in one vector M t = m 1 t ... m n t we get d M t d t = − D D O M . Δ M t (5) Where Δ = δ i , j 1 ≤ i , j ≤ n is the Laplacian matrix of the graph G E , N defined by: δ i , j = d e g i i f i = j δ i , j = − 1 i f i ≠ j ∧ i , j ∈ E δ i , j = 0 o t h e r w i s e where d e g ( i ) = ∑ ​ j : ( i , j ) ∈ E 1 is the number of adjacent nodes to the node i. Equation  is called Graph Diffusion Equation (GDE) of G N , E . Numerical schemes For temporal discretization, we employed two different schemes: Explicit scheme (forward Euler) m i ( k + 1 ) = m i ( k ) − ∑ ​ j : ( i , j ) ∈ E D D O M δ t ( m i ( k ) − m j ( k ) ) (6) Implicit scheme (backward Euler). m i ( k ) = m i ( k + 1 ) + ∑ ​ j : ( i , j ) ∈ E D D O M δ t ( m i ( k + 1 ) − m j ( k + 1 ) ) (7) The explicit scheme (Eq. ) offers simplicity but may suffer from stability issues, while the implicit scheme (Eq. ), solved using the conjugate gradient method, provides robustness at the cost of computational complexity. The detailed formulations of the numerical schemes are provided in the supplementary material S1 File. where k is the iteration number (i.e., the discrete-time index), and δ t is the time step. 2.3 Modelling microbial decomposition of organic matter processes in the graph of connected voxels Microbial activity is governed not only by the diffusion of various compounds within the soil’s pores but also by transformation processes. Dissolved Organic Matter (DOM) arises from the decomposition of both slow-decomposing Soil Organic Matter (SOM) and fast-decomposing Fresh Organic Matter (FOM). Microorganisms grow through the assimilation of DOM and breath by producing CO2. Upon death, they are transformed back into organic matter. We recall that for each voxel at the position p i : x i 1 t is the mass of MB, x i 2 t is the mass of DOM, x i 3 t is the mass of SOM, x i 4 t is the mass of FOM, x i 5 t is the mass of CO 2 . Let X j t = x 1 j t , ⋯ , x n j t represent the distribution of the j th biological parameter in the pore space. For instance, X 1 t = x 1 1 t , ⋯ , x n 1 t is the mass distribution of microorganisms in the pore space. Let v i ∈ V be a node of the graph. We model the growth of microorganisms by consuming the available dissolved organic matter according to the Monod equation. Exponential models are employed to represent both death and mass loss due to respiration through mineralized organic matter emission. The total variation of x i 1 t is described by the following equation: d x i 1 t d t = − ρ . x i 1 t − μ . x i 1 t + v D O M . x i 2 t K D O M + x i 2 t . x i 1 t (8) Where ρ is the respiration rate, μ is the mortality rate, v D O M and K D O M are respectively the maximum growth rate of MB and the constant of half saturation of DOM. In the saturated pore space, dissolved organic matter comes from dead microorganisms, the transformation of soil organic matter, and fresh organic matter. As before, a portion is allocated from DOM to microorganisms according to the Monod equation. Furthermore, in addition to diffusion processes, the total variation of x i 2 t is described by the following equation: d x i 2 t d t = − D D O M . Δ X 2 t i + β . μ . x i 1 t − v D O M . x i 2 t K D O M + x i 2 t . x i 1 t + v S O M . x i 3 t + v F O M . x i 4 t (9) Where Δ = δ i , j 1 ≤ i , j ≤ n is the Laplacian matrix of the graph G V , E , Δ X 2 t i is the i th component of Δ X 2 t , D D O M is the diffusion coefficient of DOM in water, β is the proportion of MB returning to DOM (the other fraction 1 − β returns to SOM), v F O M and v S O M the hydrolysis rate of FOM and SOM. A part of soil organic matter comes from dead microorganisms so the total variation of x i 3 t is described by the following equation: d x i 3 t d t = 1 − β μ . x i 1 t − v S O M . x i 3 t (10) The density evolution of FOM ( x i 4 t ) is presented by: d x i 4 t d t = − v F O M . x i 4 t (11) The evolution of CO2 ( x i 5 t ) is governed by microorganisms breathing, which leads to the following equations: d x i 5 t d t = ρ . x i 1 t (12) Our biological model for microbial activity in soil have been validated in , and followed the approach of the SOMKO model , and MIOR model  where microbes are subjected to three physiological processes: survival, maintenance, and growth. In this approach, the CO2 flux is computed as a linear function of biomass. By gathering all the equations, we have the following system of equations: { d x i 1 ( t ) / d t = − p . x i 1 ( t ) − μ . x i 1 ( t ) + v D O M . x i 2 ( t ) K D O M + x i 2 ( t ) . x i 1 ( t ) d x i 2 ( t ) d t = D D O M . [ Δ X 2 ( t ) ] i + β . μ . x i 1 ( t ) − v D O M . x i 2 ( t ) K D O M + x i 2 ( t ) . x i 1 ( t ) + v S O M . x i 3 ( t ) + v F O M . x i 4 ( t ) d x i 3 ( t ) d t = ( 1 − β ) μ . x i 1 ( t ) − v S O M . x i 3 ( t ) d x i 4 ( t ) d t = − v F O M . x i 4 ( t ) d x i 5 ( t ) d t = ρ . x i 1 ( t ) (13) By removing the diffusion processes of DOM from the equations, and applying an explicit scheme to the first-order derivatives of the equations we obtain the following equations for the transformation processes: { x i 1 ( t + δ t ) = x i 1 ( t ) + v D O M . x i 2 ( t ) K D O M + x i 2 ( t ) . x i 1 ( t ) . δ t − ( ρ . x i 1 ( t ) + μ . x i 1 ( t ) ) . δ t x i 2 ( t + δ t ) = x i 2 ( t + δ t ) − v D O M . x i 2 ( t ) K D O M + x i 2 ( t ) . x i 1 ( t ) . δ t + ( v S O M . x i 3 ( t ) + v F O M . x i 4 ( t ) + β . μ x i 1 ( t ) ) . δ t x i 3 ( t + δ t ) = x i 3 ( t ) + ( 1 − β ) μ . x i 1 ( t ) − v S O M . x i 3 ( t ) . δ t x i 4 ( t + δ t ) = x i 4 ( t ) − v F O M . x i 4 ( t ) . δ t x i 5 ( t + δ t ) = x i 5 ( t ) − ρ . x i 1 ( t ) . δ t (14) To obtain the voxel states after a time step δ t of transformation processes, we can apply the equations synchronously or asynchronously. Given that the transformation of the nodes is independent, we have the flexibility to perform transformations sequentially or in parallel on the graph’s nodes. In the second section, we provide a comprehensive explanation of the implementation of all discussed concepts. This includes details on graph construction from the 3D image, the conjugate gradient method, its parallelization, and the implementation of both implicit and explicit schemes. Furthermore, we provide details of implementing transformation processes, whether synchronous or asynchronous. 2.4 Implementation details Generally, simulating transformation-diffusion processes in the complex geometry captured by the 3D image involves transforming the pore space in the image into a graph of connected voxels. Subsequently, each simulation iteration applies transformation processes to the biological attributes of all graph nodes, followed by the solution of the graph diffusion equation to mimic the diffusion of the compound to be diffused in the model (DOM in this study). Transforming the 3D image into a graph of connected voxels is a straightforward and parallelizable process. It can be accomplished by enumerating the valid voxels in the image (i.e., those corresponding to the pore space). Then, looping over the voxels of the image, when we encounter a valid voxel, we identify from the six neighbouring voxels  the valid ones and construct the adjacency graph. Let t ≥ 0 be a specific time, and let X i t i ∈ N denote the distribution of the nodes at time t . The attributes of the nodes, i.e., the distribution of the pore space X i t + δ t i ∈ N after a time step δ t of microbial activity, are obtained by applying the processes in the system of Eq. synchronously (flowchart S3A in S1 File) or asynchronously (flowchart S3B in S1 File). Then, the explicit scheme (Eq. ) or implicit scheme (Eq. ) is applied in order to get the distribution after diffusion processes. In the synchronous procedure for performing microbial transformation, we calculate the contribution of biological variables to uptake, respiration, mortality, and turnover for each node. Subsequently, we update the biological variables according to the logic of the system of Eq. . In contrast, the asynchronous procedure involves sequentially updating the biological variables of the nodes due to each biological process (grow, breath, death …). That means allowing microorganisms to grow by consuming DOM, after which some will die, producing CO₂ through respiration, and transforming the dead microbial biomass into DOM and SOM. Detailed flowcharts of the algorithms for both synchronous and asynchronous transformations are described in in flowcharts S3A and S3B in S1 File. When the time steps for both transformation and diffusion are identical, the processes are performed just once per time step. However, if the time steps differ between these processes, we iterate through the process (transformation or diffusion) with the smallest time step to ensure that the dynamics of the entire transformation-diffusion system remain accurate. When using the explicit scheme, the computational process is reduced to straightforward matrix multiplications, that is relatively efficient. On the other hand, implementing the implicit scheme (Eq. ) requires solving a linear system defined by a large, sparse, symmetric, and positive definite matrix. To handle such computations effectively, we employ the preconditioned conjugate gradient method, where the conditioning matrix t i , j 0 ≤ i , j ≤ n , is defined by t i , j = 1 b i , j i f i = j , t i , j = 0 e l s e where B = b i , j 0 ≤ i , j ≤ n is the matrix of the implicit scheme defined in supplementary material S1 File. The Preconditioned Conjugate Gradient (PCG) method is a variant of the Conjugate Gradient (CG) method used for solving linear systems of equations A x = b , where A is a symmetric positive definite (SPD) matrix. The PCG method incorporates a preconditioner matrix T to improve convergence speed. The implementation was done using the C language due to its capacity for handling such complex problems: efficient memory management, computational speed required to address some aspects of the simulation like: graph construction from the 3D image, solving the graph diffusion equation through the numerical schemes, and running expensive simulations on the constructed graph. The C implementation and the image data used in this study can be found at the following links: https://github.com/mouadklai/VGA_microbial_decomposition https://www.kaggle.com/datasets/mouadklai/threed-image-sandy-loam-soil/data 2.5 Diffusional conductance coefficients approximation: theoretical framework In pore network geometrical modelling, pore space is approximated with a minimal set of maximal geometrical primitives (balls, ellipsoid, cylinders …). Then, a valuated graph is constructed from the set, where the adjacency is the geometrical adjacency between the geometrical primitives. For each pair of geometrically adjacent primitives, we need to determine the conductance, which defines the portion of mass flow to account for when simulating diffusion according to Fick’s law. In this subsection we discuss a possible way to determine approximatively the diffusional conductance coefficients of a pore network geometrical model. Let G V G A N V G A , E V G A be a voxel graph representation of a 3D image representing a pore space, where N V G A = v 1 , … , v n represents the valid voxels of the image, n is their number, and E V G A = i , j : i , j ∈ 1 , … , n ∧ v i ∩ v j ≠ ∅ encodes the adjacency between them. Let G P N G M N P N G M , E P N G M be the graph of connected geometrical primitives covering the pore space of the 3D image, where N P N G M = P 1 , … , P q is the set of geometrical primitives, q is their number, and E P N G M = i , j : i , j ∈ 1 , … , n ∧ P i ∩ P j ≠ ∅ encodes the adjacency between them. For all k ∈ 1 , … , q let V P k = v i ∈ N V G A : v i ∈ P k be the set of voxels contained within the primitive P k . For theoretical explanation, suppose that: k 1 , k 2 ∈ 1 , … , q , V P k 1 ∩ V P k 2 = ∅ k ∈ 1 , … , q V P k = N V G A Let t ≥ 0 be a time, and M P N G M t = m 1 P N G M t , … , m q P N G M t be a mass distribution of the geometrical primitives P 1 , … , P q . According to Fick’s Law of diffusion, the flow of mass between two adjacent primitives P i and P j , of volume v i and v j respectively, at time t + δ t is given by: F i , j = − D c . δ t . α i , j m i t v i − m j t v j (15) where D c is the diffusion coefficient, and α i , j is what we call the diffusional conductance between the two primitives. This conductance typically depends on factors such as the contact surface between them, the distance between their centers of mass, the form of the geometrical primitives, and the difference in volume between the two primitives. Then, the variation of mass distribution d M P N G M t d t can be derived in the same way as we have done in Section 1.2, and it is given by:Top of Form ∀ i ∈ 1 , … , q , d m i P N G M t d t = ∑ j ∈ ϑ i − D c . α i , j m i P N G M t v i − m j P N G M t v j (16) where ϑ i = j ∈ 1 , … , q : P i ∩ P j ≠ ∅ is the indexes of geometrically adjacent primitives to the primitive P i . Numerically, solving the P N G M d i f f u s i o n m o d e l , can be done by discretizing time using forward (explicit) or backward (implicit) Euler schemes. Using the explicit scheme, we get: ∀ i ∈ 1 , … , q , m i P N G M t + δ t = m i P N G M t − D c . δ t . ∑ j ∈ ϑ i . α i , j m i P N G M t v i − m j P N G M t v j (17) Then, we have ∀ i ∈ 1 , … , q , m i P N G M k + 1 = m i P N G M k − D c . δ t . ∑ j ∈ ϑ i . α i , j m i P N G M k v i − m P N G M k j v j (18) Using the implicit scheme, we get: ∀ i ∈ 1 , … , q , m i P N G M t = m i P N G M t + δ t + D c . δ t . ∑ j ∈ ϑ i . α i , j m i P N G M t + δ t v i − m j P N G M t + δ t v j (19) Then, we have ∀ i ∈ 1 , … , q , m P N G M k i = m P N G M k + 1 i + D c . δ t . ∑ j ∈ ϑ i . α i , j m P N G M k + 1 i v i − m P N G M k + 1 j v j (20) where k is the iteration number. The accuracy of the P N G M d i f f u s i o n m o d e l depends on the exactitude of calculating the coefficients α i , j : i , j ∈ E P N M . In  and  the coefficients α i , j were set to α S i , j d i , j , where S i , j represents the contact surface between the primitive P i and the adjacent primitive P j , d i , j is the distance between the centers of gravitation of the spheres i and j , and α is a unified coefficient for all adjacent primitives in the network. The value of α was determined to maximize the intercorrelation (Cosinus) c o s C L B M , C P N G M α between C L B M and C P N G M α , which are the simulated curves of mass distribution in the z-planes obtained using LBM and PNGM, respectively. c o s C L B M , C P N G M α = C L B M . C P N G M α C L B M . C P N G M α = ∑ i = 1 i = 512 C L B M i . C P N G M α i ∑ i = 1 i = 512 C L B M i 2 . ∑ i = 1 i = 512 C P N G M α i 2 (21) The maximum value of c o s C L B M , C P N G M α was achieved at α =  0.6 in the experiment described in section 3.1, using the dataset presented in this study, that was also employed in . The curve labeled “old PNGM” in represents the results obtained using PNGM. The curve obtained using LBM is not shown, as it closely resembles the one generated using VGA, which is displayed in . In the following, we present a detailed approach to approximate the coefficients using a machine learning approach and data generated using the voxel graph-based approach discussed before. Let, M V G A t = m 1 V G A t , … , m n V G A t be the voxel distribution that correspond to the primitive’s distribution M P N M t at time t, obtained by the voxel description V P k : k ∈ 1 , … , q of the primitives { P k : k ∈ 1 , … , q } . The evolution of M V G A t over time can be calculated using the graph diffusion equation through the use of one of the schemes: the explicit (equation 16) or the implicit (equation 18). The framework for the simulation of M V G A t over time was discussed in section 2.2 and validated in section 3.1. Let X p , Y p : p ∈ 1 , … , l be a dataset of couple of distributions: a random distribution X p of mass in the pore network and the corresponding next distribution Y p after a time step δ t of diffusion calculated using the voxel graph-based approach from the corresponding voxel description. Rewriting the scheme 3 for the dataset we have ∀ p ∈ 1 , … , l , ∀ i ∈ 1 , … , q , y i p = x i p − D c . δ t . ∑ j ∈ ϑ i . α i , j x i p v i − x i p v j (22) Finding α i , j : i , j ∈ E P N G M is equivalent to finding the set of parameters Θ = θ i , j : i , j ∈ E P N G M that minimize for all data X p , Y p : p ∈ 1 , … , l one of the following objectives: L 1 Θ = 1 l ∑ p = 1 l 1 q ∑ i = 1 q y i p − x i p + D c . δ t . ∑ j ∈ ϑ i . θ i , j x i p v i − x i p v j 2 (23) and L 2 Θ = 1 l ∑ p = 1 l 1 q ∑ i = 1 q x i p − y i p − D c . δ t . ∑ j ∈ ϑ i . θ i , j y i p v i − y i p v j 2 (24) The first objective calculates the L2 loss between the output of the explicit scheme (Eq.18) and the target (next time distribution) obtained using VGA. The second objective minimizes the error between the input and second member of the implicit scheme (Eq. 20) applied to the target (next time distribution) obtained using VGA. Stochastic Gradient Descent (SGD) is an efficient optimization algorithm widely used in machine learning to minimize loss functions. Unlike standard Gradient Descent, which computes gradients using the entire dataset, SGD updates parameters more frequently by using a single data point or a small batch, making it faster but introducing noise. The gradient for an objective function involving a dataset X p , Y p : p ∈ 1 , … , l (with l representing batch size) is computed as follows d d θ i , j L 1 Θ = 1 l ∑ p = 1 l 2 q D c . δ t x i p v i − x j p v j × y i p − x i p + D c . δ t ∑ f ∈ ϑ i θ i , f . x i p v i − x f p v f (25) And the gradient of the second objective is calculated using the following: d d θ i , j L 2 Θ = 1 l ∑ p = 1 l 2 q D c . δ t y i p v i − y j p v j × x i p − y i p − D c . δ t ∑ f ∈ ϑ i θ i , f . y i p v i − y f p v f (26) In each training step, we calculate the set d d θ i , j L Θ : i , j ∈ E P N M , and we update the parameters using the following stochastic gradient descent policy: ∀ i , j ∈ E P N G M , θ i , j = θ i , j − l r d d θ i , j L Θ d d θ i , j L Θ . (27) The use of the normalized gradient d d θ i , j L Θ d d θ i , j L Θ instead of the gradient d d θ i , j L Θ addresses issues related to extremely small gradient magnitudes. This normalization helps in maintaining numerical stability, ensuring consistent and meaningful updates to the parameters, and potentially improving the convergence rate of the optimization process. Using the attained value of the objective function L Θ at each training iteration, we decide whether to stop training, continue, or adjust the learning rate ( l r ) and restart training. After the training using stochastic gradient descent, we extract the learned parameters that will approximate the coefficient α i , j : i , j ∈ E P N G M . Note that the estimated diffusional conductance coefficients are theoretically applicable across all diffusion scenarios, regardless of the chosen diffusion coefficient or time step used in training, as they remain independent of these parameters. In results and discussion section we test the discussed theoretical framework in predicting the diffusional conductance coefficient for the ball network model discussed before. 2.6 Dataset description To compare the performance of the implemented numerical schemes, we used 3D microcomputed tomography (µCT) image data obtained from an experimental sandy loam soil collected at the Bullion Field site, located at the James Hutton Institute in Invergowrie, Scotland. The dataset was previously used in studies such as  ]. 2.6.1 Soil and imaging characteristics. The Bullion Field soil sample is classified as sandy loam, comprising 71% sand, 19% silt, and 10% clay by mass. This soil texture is ideal for investigating microbial activity due to its heterogeneous pore structure. The µ CT scan was performed at a pixel resolution of 24 µm, which is sufficiently high to capture the range of pore sizes associated with microbial habitats and nutrient diffusion. A cubic sub-volume of 512 3 voxels was extracted from the 3D stack, corresponding to a physical volume of approximately 1.855 cm³. 2.6.2 Image segmentation and pore structure analysis. To differentiate the solid and pore phases within the soil, we applied the indicator kriging method , yielding a segmented 3D binary image. The segmentation reveals visible 17% and 8% porosities for soil bulk densities of 1.2 g cm ⁻ ³ and 1.6 g cm ⁻ ³, respectively. The segmented image was analyzed to extract pore network characteristics, particularly focusing on maximal inscribed spheres within the pore space, following the methodology outlined in . This approach ensures the extracted pore network accurately represents the critical pathways for microbial movement and nutrient diffusion. 2.6.3 Boundary conditions and simulation setup. In the simulations, the soil sample was treated as a closed environment, which is a critical consideration for the boundary conditions applied. This assumption was made to reflect the controlled conditions typically found in experimental settings, allowing for more precise comparisons of microbial dynamics and diffusion processes within the soil structure. and in the manuscript illustrate representative cross-sections and 3D views of the binary image, highlighting the heterogeneous pore structure. For pore network geometrical modeling, depicts a view of the ball network wherein we specifically choose the balls whose centers reside within the region bounded by 50 , 100 , 50 , 100 , 150 , 200 . 

In this study, we only considered the condition that soil is fully saturated, i.e., soil pore space is filled with water. We model the dynamics of microbial decomposition using five key compounds, as illustrated in : Microbial biomass (MB): Represents the mass of microorganisms in the sample. Microbial respiration (CO2): Indicates the produced carbon dioxide through microbial decomposition, indicative of microbial growth. 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 the hydrolysis of FOM and SOM and biomass recycling, available for microbial uptake or transport within the pore space. Note again that for environments where saturation is not complete, the same simulation framework can be used after a method of drainage  which refers to the method used to redistribute liquid and air within a porous medium, typically by gradually removing liquid from the system. This process mimics the movement of liquid due to capillary forces and gravity, leading to the formation of air-filled pores. In simulations, this is often achieved by progressively reducing the liquid content and adjusting the pore space accordingly. This process allows the model to represent partially saturated conditions before transforming the system into a two-phase configuration by treating the air phase as a non-interacting solid phase . Let I i , j , k i , j , k be a 3D binary image where the voxels forming pore space (void voxels) are tagged by 0 and the voxels attached to solid matter (solid voxels) are tagged by 1. Let V = i , j , k : I i , j , k = 0 = v 1 , … , v n be the set of pore space voxels and N = 1 , … , n the index set of V , where n = c a r d V represents the number of voxels of the pore space. We construct an adjacency valuated graph G V , E where V is the set of nodes, and E = { i , j ∈ N 2 | v i ∩ v j ≠ ∅ } is the set of edges. In this work, we use 6-connecvity between voxels . Let t be the variable representing time. To each voxel v i of the pore space, we attach: the mass of MB ( x i 1 t , the mass of DOM ( x i 2 t , the mass of SOM x || i 3 t , the mass of FOM ( x i 4 t , the mass of CO 2 x || i 5 t . To each node v i let us denote by X i t = x i 1 t , x i 2 t , x i 3 t , x i 4 t , x i 5 t the masses of the compounds contained within voxel v i at time t. We simulate microbial activity by updating the constructed graph according to transformation and diffusion processes described in . The principle is to break down the complex processes, that can be mathematically modelled, into simpler steps. Let us assume that we get a representation of the pore space and its compounds at time t in the form of G t V , E , X i t i ∈ N . To get G t + δ t V , E , X i t + δ t i ∈ N at time t + δ t , after microbial activity processes we apply a transformation model that encodes the biological conversion of different compounds within each voxel. Next, we employ a discretized model to mimic the diffusion of compounds between the connected voxels. The modeling approach is based on dynamics splitting approach, where the diffusion and reaction processes are handled separately. The spatial discretization is fixed according to the computed tomography image resolution. The same framework can be used to model diffusion-transformation processes at any resolution, as long as the parameters are converted based on the voxel size, which depends on the resolution of the 3D image. 

2.2.1 Fick’s law of diffusion. Fick’s law of diffusion describes the process of molecular or particle movement through a medium, such as gases, liquids, or solids. It is expressed as follows : J = − D d C d x (1) Where J represents the diffusion flux (amount of substance per unit area per unit time), D is the of the substance and d c / d x the concentration gradient along the diffusion direction. 2.2.2 Diffusion modelling: graph diffusion equation. We consider the diffusion of dissolved organic matter (DOM). Let m i t = x i 2 t represent the mass of DOM at time t in the voxel p i . The diffusion process in our framework is modeled using a graph-based approach where voxels represent graph nodes, and edges connect neighboring voxels. The concentration at a voxel corresponds to its mass since the voxel volume is normalized. The mass flow between time t and time t + δ t from the voxel at position p i to its neighboring voxel at position p i is given by the first Fick’s law as follows: J i , j = − D D O M . m i t − m j t . δ t (2) where D D O M is the DOM diffusion coefficient in water. Thus, the DOM mass at voxel i at time t + δ t due to mass exchange with neighbouring voxels is expressed as follows: m i t + δ t − m i t = ∑ j : i , j ∈ E − D D O M δ m i , j t δ t (3) Then the governing equation is: d m i t d t = ∑ j : i , j ∈ E − D D O M m i t − m j t (4) Gathering all the masses of all the voxels of the pore space in one vector M t = m 1 t ... m n t we get d M t d t = − D D O M . Δ M t (5) Where Δ = δ i , j 1 ≤ i , j ≤ n is the Laplacian matrix of the graph G E , N defined by: δ i , j = d e g i i f i = j δ i , j = − 1 i f i ≠ j ∧ i , j ∈ E δ i , j = 0 o t h e r w i s e where d e g ( i ) = ∑ ​ j : ( i , j ) ∈ E 1 is the number of adjacent nodes to the node i. Equation  is called Graph Diffusion Equation (GDE) of G N , E . Numerical schemes For temporal discretization, we employed two different schemes: Explicit scheme (forward Euler) m i ( k + 1 ) = m i ( k ) − ∑ ​ j : ( i , j ) ∈ E D D O M δ t ( m i ( k ) − m j ( k ) ) (6) Implicit scheme (backward Euler). m i ( k ) = m i ( k + 1 ) + ∑ ​ j : ( i , j ) ∈ E D D O M δ t ( m i ( k + 1 ) − m j ( k + 1 ) ) (7) The explicit scheme (Eq. ) offers simplicity but may suffer from stability issues, while the implicit scheme (Eq. ), solved using the conjugate gradient method, provides robustness at the cost of computational complexity. The detailed formulations of the numerical schemes are provided in the supplementary material S1 File. where k is the iteration number (i.e., the discrete-time index), and δ t is the time step. 

Fick’s law of diffusion describes the process of molecular or particle movement through a medium, such as gases, liquids, or solids. It is expressed as follows : J = − D d C d x (1) Where J represents the diffusion flux (amount of substance per unit area per unit time), D is the of the substance and d c / d x the concentration gradient along the diffusion direction. 

We consider the diffusion of dissolved organic matter (DOM). Let m i t = x i 2 t represent the mass of DOM at time t in the voxel p i . The diffusion process in our framework is modeled using a graph-based approach where voxels represent graph nodes, and edges connect neighboring voxels. The concentration at a voxel corresponds to its mass since the voxel volume is normalized. The mass flow between time t and time t + δ t from the voxel at position p i to its neighboring voxel at position p i is given by the first Fick’s law as follows: J i , j = − D D O M . m i t − m j t . δ t (2) where D D O M is the DOM diffusion coefficient in water. Thus, the DOM mass at voxel i at time t + δ t due to mass exchange with neighbouring voxels is expressed as follows: m i t + δ t − m i t = ∑ j : i , j ∈ E − D D O M δ m i , j t δ t (3) Then the governing equation is: d m i t d t = ∑ j : i , j ∈ E − D D O M m i t − m j t (4) Gathering all the masses of all the voxels of the pore space in one vector M t = m 1 t ... m n t we get d M t d t = − D D O M . Δ M t (5) Where Δ = δ i , j 1 ≤ i , j ≤ n is the Laplacian matrix of the graph G E , N defined by: δ i , j = d e g i i f i = j δ i , j = − 1 i f i ≠ j ∧ i , j ∈ E δ i , j = 0 o t h e r w i s e where d e g ( i ) = ∑ ​ j : ( i , j ) ∈ E 1 is the number of adjacent nodes to the node i. Equation  is called Graph Diffusion Equation (GDE) of G N , E . Numerical schemes For temporal discretization, we employed two different schemes: Explicit scheme (forward Euler) m i ( k + 1 ) = m i ( k ) − ∑ ​ j : ( i , j ) ∈ E D D O M δ t ( m i ( k ) − m j ( k ) ) (6) Implicit scheme (backward Euler). m i ( k ) = m i ( k + 1 ) + ∑ ​ j : ( i , j ) ∈ E D D O M δ t ( m i ( k + 1 ) − m j ( k + 1 ) ) (7) The explicit scheme (Eq. ) offers simplicity but may suffer from stability issues, while the implicit scheme (Eq. ), solved using the conjugate gradient method, provides robustness at the cost of computational complexity. The detailed formulations of the numerical schemes are provided in the supplementary material S1 File. where k is the iteration number (i.e., the discrete-time index), and δ t is the time step. 

Microbial activity is governed not only by the diffusion of various compounds within the soil’s pores but also by transformation processes. Dissolved Organic Matter (DOM) arises from the decomposition of both slow-decomposing Soil Organic Matter (SOM) and fast-decomposing Fresh Organic Matter (FOM). Microorganisms grow through the assimilation of DOM and breath by producing CO2. Upon death, they are transformed back into organic matter. We recall that for each voxel at the position p i : x i 1 t is the mass of MB, x i 2 t is the mass of DOM, x i 3 t is the mass of SOM, x i 4 t is the mass of FOM, x i 5 t is the mass of CO 2 . Let X j t = x 1 j t , ⋯ , x n j t represent the distribution of the j th biological parameter in the pore space. For instance, X 1 t = x 1 1 t , ⋯ , x n 1 t is the mass distribution of microorganisms in the pore space. Let v i ∈ V be a node of the graph. We model the growth of microorganisms by consuming the available dissolved organic matter according to the Monod equation. Exponential models are employed to represent both death and mass loss due to respiration through mineralized organic matter emission. The total variation of x i 1 t is described by the following equation: d x i 1 t d t = − ρ . x i 1 t − μ . x i 1 t + v D O M . x i 2 t K D O M + x i 2 t . x i 1 t (8) Where ρ is the respiration rate, μ is the mortality rate, v D O M and K D O M are respectively the maximum growth rate of MB and the constant of half saturation of DOM. In the saturated pore space, dissolved organic matter comes from dead microorganisms, the transformation of soil organic matter, and fresh organic matter. As before, a portion is allocated from DOM to microorganisms according to the Monod equation. Furthermore, in addition to diffusion processes, the total variation of x i 2 t is described by the following equation: d x i 2 t d t = − D D O M . Δ X 2 t i + β . μ . x i 1 t − v D O M . x i 2 t K D O M + x i 2 t . x i 1 t + v S O M . x i 3 t + v F O M . x i 4 t (9) Where Δ = δ i , j 1 ≤ i , j ≤ n is the Laplacian matrix of the graph G V , E , Δ X 2 t i is the i th component of Δ X 2 t , D D O M is the diffusion coefficient of DOM in water, β is the proportion of MB returning to DOM (the other fraction 1 − β returns to SOM), v F O M and v S O M the hydrolysis rate of FOM and SOM. A part of soil organic matter comes from dead microorganisms so the total variation of x i 3 t is described by the following equation: d x i 3 t d t = 1 − β μ . x i 1 t − v S O M . x i 3 t (10) The density evolution of FOM ( x i 4 t ) is presented by: d x i 4 t d t = − v F O M . x i 4 t (11) The evolution of CO2 ( x i 5 t ) is governed by microorganisms breathing, which leads to the following equations: d x i 5 t d t = ρ . x i 1 t (12) Our biological model for microbial activity in soil have been validated in , and followed the approach of the SOMKO model , and MIOR model  where microbes are subjected to three physiological processes: survival, maintenance, and growth. In this approach, the CO2 flux is computed as a linear function of biomass. By gathering all the equations, we have the following system of equations: { d x i 1 ( t ) / d t = − p . x i 1 ( t ) − μ . x i 1 ( t ) + v D O M . x i 2 ( t ) K D O M + x i 2 ( t ) . x i 1 ( t ) d x i 2 ( t ) d t = D D O M . [ Δ X 2 ( t ) ] i + β . μ . x i 1 ( t ) − v D O M . x i 2 ( t ) K D O M + x i 2 ( t ) . x i 1 ( t ) + v S O M . x i 3 ( t ) + v F O M . x i 4 ( t ) d x i 3 ( t ) d t = ( 1 − β ) μ . x i 1 ( t ) − v S O M . x i 3 ( t ) d x i 4 ( t ) d t = − v F O M . x i 4 ( t ) d x i 5 ( t ) d t = ρ . x i 1 ( t ) (13) By removing the diffusion processes of DOM from the equations, and applying an explicit scheme to the first-order derivatives of the equations we obtain the following equations for the transformation processes: { x i 1 ( t + δ t ) = x i 1 ( t ) + v D O M . x i 2 ( t ) K D O M + x i 2 ( t ) . x i 1 ( t ) . δ t − ( ρ . x i 1 ( t ) + μ . x i 1 ( t ) ) . δ t x i 2 ( t + δ t ) = x i 2 ( t + δ t ) − v D O M . x i 2 ( t ) K D O M + x i 2 ( t ) . x i 1 ( t ) . δ t + ( v S O M . x i 3 ( t ) + v F O M . x i 4 ( t ) + β . μ x i 1 ( t ) ) . δ t x i 3 ( t + δ t ) = x i 3 ( t ) + ( 1 − β ) μ . x i 1 ( t ) − v S O M . x i 3 ( t ) . δ t x i 4 ( t + δ t ) = x i 4 ( t ) − v F O M . x i 4 ( t ) . δ t x i 5 ( t + δ t ) = x i 5 ( t ) − ρ . x i 1 ( t ) . δ t (14) To obtain the voxel states after a time step δ t of transformation processes, we can apply the equations synchronously or asynchronously. Given that the transformation of the nodes is independent, we have the flexibility to perform transformations sequentially or in parallel on the graph’s nodes. In the second section, we provide a comprehensive explanation of the implementation of all discussed concepts. This includes details on graph construction from the 3D image, the conjugate gradient method, its parallelization, and the implementation of both implicit and explicit schemes. Furthermore, we provide details of implementing transformation processes, whether synchronous or asynchronous. 

Generally, simulating transformation-diffusion processes in the complex geometry captured by the 3D image involves transforming the pore space in the image into a graph of connected voxels. Subsequently, each simulation iteration applies transformation processes to the biological attributes of all graph nodes, followed by the solution of the graph diffusion equation to mimic the diffusion of the compound to be diffused in the model (DOM in this study). Transforming the 3D image into a graph of connected voxels is a straightforward and parallelizable process. It can be accomplished by enumerating the valid voxels in the image (i.e., those corresponding to the pore space). Then, looping over the voxels of the image, when we encounter a valid voxel, we identify from the six neighbouring voxels  the valid ones and construct the adjacency graph. Let t ≥ 0 be a specific time, and let X i t i ∈ N denote the distribution of the nodes at time t . The attributes of the nodes, i.e., the distribution of the pore space X i t + δ t i ∈ N after a time step δ t of microbial activity, are obtained by applying the processes in the system of Eq. synchronously (flowchart S3A in S1 File) or asynchronously (flowchart S3B in S1 File). Then, the explicit scheme (Eq. ) or implicit scheme (Eq. ) is applied in order to get the distribution after diffusion processes. In the synchronous procedure for performing microbial transformation, we calculate the contribution of biological variables to uptake, respiration, mortality, and turnover for each node. Subsequently, we update the biological variables according to the logic of the system of Eq. . In contrast, the asynchronous procedure involves sequentially updating the biological variables of the nodes due to each biological process (grow, breath, death …). That means allowing microorganisms to grow by consuming DOM, after which some will die, producing CO₂ through respiration, and transforming the dead microbial biomass into DOM and SOM. Detailed flowcharts of the algorithms for both synchronous and asynchronous transformations are described in in flowcharts S3A and S3B in S1 File. When the time steps for both transformation and diffusion are identical, the processes are performed just once per time step. However, if the time steps differ between these processes, we iterate through the process (transformation or diffusion) with the smallest time step to ensure that the dynamics of the entire transformation-diffusion system remain accurate. When using the explicit scheme, the computational process is reduced to straightforward matrix multiplications, that is relatively efficient. On the other hand, implementing the implicit scheme (Eq. ) requires solving a linear system defined by a large, sparse, symmetric, and positive definite matrix. To handle such computations effectively, we employ the preconditioned conjugate gradient method, where the conditioning matrix t i , j 0 ≤ i , j ≤ n , is defined by t i , j = 1 b i , j i f i = j , t i , j = 0 e l s e where B = b i , j 0 ≤ i , j ≤ n is the matrix of the implicit scheme defined in supplementary material S1 File. The Preconditioned Conjugate Gradient (PCG) method is a variant of the Conjugate Gradient (CG) method used for solving linear systems of equations A x = b , where A is a symmetric positive definite (SPD) matrix. The PCG method incorporates a preconditioner matrix T to improve convergence speed. The implementation was done using the C language due to its capacity for handling such complex problems: efficient memory management, computational speed required to address some aspects of the simulation like: graph construction from the 3D image, solving the graph diffusion equation through the numerical schemes, and running expensive simulations on the constructed graph. The C implementation and the image data used in this study can be found at the following links: https://github.com/mouadklai/VGA_microbial_decomposition https://www.kaggle.com/datasets/mouadklai/threed-image-sandy-loam-soil/data 

In pore network geometrical modelling, pore space is approximated with a minimal set of maximal geometrical primitives (balls, ellipsoid, cylinders …). Then, a valuated graph is constructed from the set, where the adjacency is the geometrical adjacency between the geometrical primitives. For each pair of geometrically adjacent primitives, we need to determine the conductance, which defines the portion of mass flow to account for when simulating diffusion according to Fick’s law. In this subsection we discuss a possible way to determine approximatively the diffusional conductance coefficients of a pore network geometrical model. Let G V G A N V G A , E V G A be a voxel graph representation of a 3D image representing a pore space, where N V G A = v 1 , … , v n represents the valid voxels of the image, n is their number, and E V G A = i , j : i , j ∈ 1 , … , n ∧ v i ∩ v j ≠ ∅ encodes the adjacency between them. Let G P N G M N P N G M , E P N G M be the graph of connected geometrical primitives covering the pore space of the 3D image, where N P N G M = P 1 , … , P q is the set of geometrical primitives, q is their number, and E P N G M = i , j : i , j ∈ 1 , … , n ∧ P i ∩ P j ≠ ∅ encodes the adjacency between them. For all k ∈ 1 , … , q let V P k = v i ∈ N V G A : v i ∈ P k be the set of voxels contained within the primitive P k . For theoretical explanation, suppose that: k 1 , k 2 ∈ 1 , … , q , V P k 1 ∩ V P k 2 = ∅ k ∈ 1 , … , q V P k = N V G A Let t ≥ 0 be a time, and M P N G M t = m 1 P N G M t , … , m q P N G M t be a mass distribution of the geometrical primitives P 1 , … , P q . According to Fick’s Law of diffusion, the flow of mass between two adjacent primitives P i and P j , of volume v i and v j respectively, at time t + δ t is given by: F i , j = − D c . δ t . α i , j m i t v i − m j t v j (15) where D c is the diffusion coefficient, and α i , j is what we call the diffusional conductance between the two primitives. This conductance typically depends on factors such as the contact surface between them, the distance between their centers of mass, the form of the geometrical primitives, and the difference in volume between the two primitives. Then, the variation of mass distribution d M P N G M t d t can be derived in the same way as we have done in Section 1.2, and it is given by:Top of Form ∀ i ∈ 1 , … , q , d m i P N G M t d t = ∑ j ∈ ϑ i − D c . α i , j m i P N G M t v i − m j P N G M t v j (16) where ϑ i = j ∈ 1 , … , q : P i ∩ P j ≠ ∅ is the indexes of geometrically adjacent primitives to the primitive P i . Numerically, solving the P N G M d i f f u s i o n m o d e l , can be done by discretizing time using forward (explicit) or backward (implicit) Euler schemes. Using the explicit scheme, we get: ∀ i ∈ 1 , … , q , m i P N G M t + δ t = m i P N G M t − D c . δ t . ∑ j ∈ ϑ i . α i , j m i P N G M t v i − m j P N G M t v j (17) Then, we have ∀ i ∈ 1 , … , q , m i P N G M k + 1 = m i P N G M k − D c . δ t . ∑ j ∈ ϑ i . α i , j m i P N G M k v i − m P N G M k j v j (18) Using the implicit scheme, we get: ∀ i ∈ 1 , … , q , m i P N G M t = m i P N G M t + δ t + D c . δ t . ∑ j ∈ ϑ i . α i , j m i P N G M t + δ t v i − m j P N G M t + δ t v j (19) Then, we have ∀ i ∈ 1 , … , q , m P N G M k i = m P N G M k + 1 i + D c . δ t . ∑ j ∈ ϑ i . α i , j m P N G M k + 1 i v i − m P N G M k + 1 j v j (20) where k is the iteration number. The accuracy of the P N G M d i f f u s i o n m o d e l depends on the exactitude of calculating the coefficients α i , j : i , j ∈ E P N M . In  and  the coefficients α i , j were set to α S i , j d i , j , where S i , j represents the contact surface between the primitive P i and the adjacent primitive P j , d i , j is the distance between the centers of gravitation of the spheres i and j , and α is a unified coefficient for all adjacent primitives in the network. The value of α was determined to maximize the intercorrelation (Cosinus) c o s C L B M , C P N G M α between C L B M and C P N G M α , which are the simulated curves of mass distribution in the z-planes obtained using LBM and PNGM, respectively. c o s C L B M , C P N G M α = C L B M . C P N G M α C L B M . C P N G M α = ∑ i = 1 i = 512 C L B M i . C P N G M α i ∑ i = 1 i = 512 C L B M i 2 . ∑ i = 1 i = 512 C P N G M α i 2 (21) The maximum value of c o s C L B M , C P N G M α was achieved at α =  0.6 in the experiment described in section 3.1, using the dataset presented in this study, that was also employed in . The curve labeled “old PNGM” in represents the results obtained using PNGM. The curve obtained using LBM is not shown, as it closely resembles the one generated using VGA, which is displayed in . In the following, we present a detailed approach to approximate the coefficients using a machine learning approach and data generated using the voxel graph-based approach discussed before. Let, M V G A t = m 1 V G A t , … , m n V G A t be the voxel distribution that correspond to the primitive’s distribution M P N M t at time t, obtained by the voxel description V P k : k ∈ 1 , … , q of the primitives { P k : k ∈ 1 , … , q } . The evolution of M V G A t over time can be calculated using the graph diffusion equation through the use of one of the schemes: the explicit (equation 16) or the implicit (equation 18). The framework for the simulation of M V G A t over time was discussed in section 2.2 and validated in section 3.1. Let X p , Y p : p ∈ 1 , … , l be a dataset of couple of distributions: a random distribution X p of mass in the pore network and the corresponding next distribution Y p after a time step δ t of diffusion calculated using the voxel graph-based approach from the corresponding voxel description. Rewriting the scheme 3 for the dataset we have ∀ p ∈ 1 , … , l , ∀ i ∈ 1 , … , q , y i p = x i p − D c . δ t . ∑ j ∈ ϑ i . α i , j x i p v i − x i p v j (22) Finding α i , j : i , j ∈ E P N G M is equivalent to finding the set of parameters Θ = θ i , j : i , j ∈ E P N G M that minimize for all data X p , Y p : p ∈ 1 , … , l one of the following objectives: L 1 Θ = 1 l ∑ p = 1 l 1 q ∑ i = 1 q y i p − x i p + D c . δ t . ∑ j ∈ ϑ i . θ i , j x i p v i − x i p v j 2 (23) and L 2 Θ = 1 l ∑ p = 1 l 1 q ∑ i = 1 q x i p − y i p − D c . δ t . ∑ j ∈ ϑ i . θ i , j y i p v i − y i p v j 2 (24) The first objective calculates the L2 loss between the output of the explicit scheme (Eq.18) and the target (next time distribution) obtained using VGA. The second objective minimizes the error between the input and second member of the implicit scheme (Eq. 20) applied to the target (next time distribution) obtained using VGA. Stochastic Gradient Descent (SGD) is an efficient optimization algorithm widely used in machine learning to minimize loss functions. Unlike standard Gradient Descent, which computes gradients using the entire dataset, SGD updates parameters more frequently by using a single data point or a small batch, making it faster but introducing noise. The gradient for an objective function involving a dataset X p , Y p : p ∈ 1 , … , l (with l representing batch size) is computed as follows d d θ i , j L 1 Θ = 1 l ∑ p = 1 l 2 q D c . δ t x i p v i − x j p v j × y i p − x i p + D c . δ t ∑ f ∈ ϑ i θ i , f . x i p v i − x f p v f (25) And the gradient of the second objective is calculated using the following: d d θ i , j L 2 Θ = 1 l ∑ p = 1 l 2 q D c . δ t y i p v i − y j p v j × x i p − y i p − D c . δ t ∑ f ∈ ϑ i θ i , f . y i p v i − y f p v f (26) In each training step, we calculate the set d d θ i , j L Θ : i , j ∈ E P N M , and we update the parameters using the following stochastic gradient descent policy: ∀ i , j ∈ E P N G M , θ i , j = θ i , j − l r d d θ i , j L Θ d d θ i , j L Θ . (27) The use of the normalized gradient d d θ i , j L Θ d d θ i , j L Θ instead of the gradient d d θ i , j L Θ addresses issues related to extremely small gradient magnitudes. This normalization helps in maintaining numerical stability, ensuring consistent and meaningful updates to the parameters, and potentially improving the convergence rate of the optimization process. Using the attained value of the objective function L Θ at each training iteration, we decide whether to stop training, continue, or adjust the learning rate ( l r ) and restart training. After the training using stochastic gradient descent, we extract the learned parameters that will approximate the coefficient α i , j : i , j ∈ E P N G M . Note that the estimated diffusional conductance coefficients are theoretically applicable across all diffusion scenarios, regardless of the chosen diffusion coefficient or time step used in training, as they remain independent of these parameters. In results and discussion section we test the discussed theoretical framework in predicting the diffusional conductance coefficient for the ball network model discussed before. 

To compare the performance of the implemented numerical schemes, we used 3D microcomputed tomography (µCT) image data obtained from an experimental sandy loam soil collected at the Bullion Field site, located at the James Hutton Institute in Invergowrie, Scotland. The dataset was previously used in studies such as  ]. 2.6.1 Soil and imaging characteristics. The Bullion Field soil sample is classified as sandy loam, comprising 71% sand, 19% silt, and 10% clay by mass. This soil texture is ideal for investigating microbial activity due to its heterogeneous pore structure. The µ CT scan was performed at a pixel resolution of 24 µm, which is sufficiently high to capture the range of pore sizes associated with microbial habitats and nutrient diffusion. A cubic sub-volume of 512 3 voxels was extracted from the 3D stack, corresponding to a physical volume of approximately 1.855 cm³. 2.6.2 Image segmentation and pore structure analysis. To differentiate the solid and pore phases within the soil, we applied the indicator kriging method , yielding a segmented 3D binary image. The segmentation reveals visible 17% and 8% porosities for soil bulk densities of 1.2 g cm ⁻ ³ and 1.6 g cm ⁻ ³, respectively. The segmented image was analyzed to extract pore network characteristics, particularly focusing on maximal inscribed spheres within the pore space, following the methodology outlined in . This approach ensures the extracted pore network accurately represents the critical pathways for microbial movement and nutrient diffusion. 2.6.3 Boundary conditions and simulation setup. In the simulations, the soil sample was treated as a closed environment, which is a critical consideration for the boundary conditions applied. This assumption was made to reflect the controlled conditions typically found in experimental settings, allowing for more precise comparisons of microbial dynamics and diffusion processes within the soil structure. and in the manuscript illustrate representative cross-sections and 3D views of the binary image, highlighting the heterogeneous pore structure. For pore network geometrical modeling, depicts a view of the ball network wherein we specifically choose the balls whose centers reside within the region bounded by 50 , 100 , 50 , 100 , 150 , 200 . 

The Bullion Field soil sample is classified as sandy loam, comprising 71% sand, 19% silt, and 10% clay by mass. This soil texture is ideal for investigating microbial activity due to its heterogeneous pore structure. The µ CT scan was performed at a pixel resolution of 24 µm, which is sufficiently high to capture the range of pore sizes associated with microbial habitats and nutrient diffusion. A cubic sub-volume of 512 3 voxels was extracted from the 3D stack, corresponding to a physical volume of approximately 1.855 cm³. 

To differentiate the solid and pore phases within the soil, we applied the indicator kriging method , yielding a segmented 3D binary image. The segmentation reveals visible 17% and 8% porosities for soil bulk densities of 1.2 g cm ⁻ ³ and 1.6 g cm ⁻ ³, respectively. The segmented image was analyzed to extract pore network characteristics, particularly focusing on maximal inscribed spheres within the pore space, following the methodology outlined in . This approach ensures the extracted pore network accurately represents the critical pathways for microbial movement and nutrient diffusion. 

In the simulations, the soil sample was treated as a closed environment, which is a critical consideration for the boundary conditions applied. This assumption was made to reflect the controlled conditions typically found in experimental settings, allowing for more precise comparisons of microbial dynamics and diffusion processes within the soil structure. and in the manuscript illustrate representative cross-sections and 3D views of the binary image, highlighting the heterogeneous pore structure. For pore network geometrical modeling, depicts a view of the ball network wherein we specifically choose the balls whose centers reside within the region bounded by 50 , 100 , 50 , 100 , 150 , 200 . 

In this section, we compare simulations using the presented approach that we call the Voxel Graph-based Approach (VGA), with LBM-based simulation and pore network-based simulations. For LBM-based simulation, we use the approach outlined in , where diffusion is computed from the lattice-Boltzmann equation, and transformation is calculated using the synchronous scheme. For the pore network geometrical model, we use the generic tool discussed in . This tool models the complex shape of the pore space using a graph of connected balls, and formulate diffusion using the first Fick’s law. The mass transport between the balls was calibrated to fit the maximum intercorrelation between simulated mass distributions after a simulation time using the ball approach and LBM approach . 3.1 Comparison of diffusion simulations: VGA, LBM and PNGM-based simulations To compare diffusion simulation using the graph diffusion equation with simulations using the Lattice Boltzmann method and the pore network model, we draw upon a prior experience detailed in , which was conducted to compare diffusion using a pore network model (ball network model) and LBM, and we outline our approach as follows: Uniformly, we distributed M 0 = 592.7593 m g of DOM on the first two z layers ( z = 1 and z = 2 ) of the 512 × 512 × 512 image. In our approach and the LBM, we selected all the voxels of the two layers ( z = 0 and z = 1 ) corresponding to the pore space and distributed M 0 uniformly among the selected voxels. For the ball network model, we selected the balls intersecting the two layers and distributed the mass M 0 homogeneously across them (i.e., applying the same concentration to all the balls). The experiment aims to compare mass profiles in each layer of the image after some time of 1.76 h o u r s . We ran the simulations using the three methods and obtained the final distributions, from which we calculated the mass in each layer of the image. Subsequently, we plot these results together for comparison. The basic principle of pore network geometrical modelling (PNGM) is to globally transport mass from one ball to the connected balls following Fick’s laws. Therefore, we need to calibrate the diffusional conductance. We employed a general coefficient α = 0.6 , which we multiplied by the contact surface area between all pair of geometrically connected balls, calibrated from LBM . Details on how the coefficient α is determined are provided in Section 2.4. illustrates a global comparison of the overall accuracy of each method. PNGM was simulated using the implicit scheme with a time step of 15 seconds, while GDE was solved using the implicit scheme with a time step of 30 seconds. In contrast, the time step for LBM is 0.43 seconds, calculated based on the 3D image resolution. The results obtained with the implicit scheme of the GDE are very close to those of LBM and much more accurate than those of the implicit scheme of PNGM. The use of voxels to discretize the space is similar to LBM; however, our approach models diffusion as mass transport according to Fick’s law, whereas LBM uses particle-based collision and propagation processes. Although both methods are conceptually similar, our approach involves less calculations. In contrast, PNGM uses geometrical primitives that are larger than a single voxel and vary in size. This approach makes modeling mass transport in the pore space using Fick’s law less precise. The explicit scheme for solving the graph diffusion equation requires the use of small-time steps to maintain stability and accuracy in the numerical solution. This necessity arises because each time step needs to be small enough to ensure that the diffusion process is accurately captured and does not lead to instability or inaccuracies in the results. As a result, the explicit scheme involves frequent updates and calculations, which can become computationally expensive, especially for large-scale problems or long simulation times. In contrast, the implicit scheme offers a more computationally efficient alternative by allowing the use of extended time steps. This approach improves computational efficiency because it reduces the number of required updates and calculations, as larger time steps cover a greater extent of the diffusion process in each iteration. Although implicit methods are typically more complex to implement, they often provide better stability and accuracy over longer time periods. illustrates the comparison between these two schemes. In the explicit scheme, a time step of 0.1 seconds is used, resulting in frequent updates and higher computational costs. On the other hand, the implicit scheme employs a relatively extended time step of 30 seconds, demonstrating a reduction in computational effort while still maintaining accuracy in the simulation results. 3.2 Comparison of microbial decomposition simulations: VGA, LBM, and PNGM-based simulations In this section, we validate the efficiency of our method for simulating microbial decomposition of organic matter by comparing its results with those obtained using the LBM-based approach and the PNGM-based method. We conducted numerical simulations of microbial decomposition using the LBM-based approach, the Voxel Graph-based Approach (VGA), and the ball network-based approach on the same dataset previously discussed. The simulation procedure is as follows: Initially, we evenly distributed M 2 0 = 289.5 μ g of DOM and introduced 1000 bacterial spots (clusters) into the pore space at random locations. This patchiness of bacterial distribution, as a result of cell growth mechanisms and environmental constraints, has been investigated in . The clusters collectively represent 5.2 × 10 7 bacteria cells, which correspond to M 1 0 = 5.2 × 10 7 × 5.41 × 10 − 8 μ g = 2.8132 μ g of carbon. For the PNGM-based approach, we distributed the mass M 2 0 homogeneously, (i.e., the same concentration in each ball), and we put bacterial mass into the balls corresponding to the random location of the spots. The total mass in the pore space initially is as follows: M 1 0 , M 2 0 , M 3 0 , M 4 0 , M 5 0 , where M 3 0 = M 4 0 = M 5 0 = 0 , corresponding to 99.0376 % of dissolved organic matter and 0.962 % of living microorganisms. The biological parameters employed in this study were adopted from  for Arthrobacter sp. 9R that was calibrated from real experimental data. These parameters are as follows: ρ = 0.2 d a y − 1 is the relative respiration rate, μ = 0.5 d a y − 1 is the relative mortality rate, β = 0.55 is the proportion of microbial mass (MB) returning to DOM (the remaining portion is for MB returning to SOM), v S O M = 0.01 d a y − 1 and v F O M = 0.3 d a y − 1 are the relative decomposition rates for SOM and FOM respectively, v D O M = 9.6 d a y − 1 and K D O M = 0.001 g C . g − 1 are the maximum relative growth rate of MB, and the half saturation constant of DOM respectively. D D O M = 100950 v o x e l 2 . d a y − 1 is the diffusion coefficient of DOM in water. The value used corresponds to the diffusion coefficient of an organic molecule in water . Parameterization of our model has been performed on experimental measurements of growth and respiration curves of Arthrobacter sp cultivated in optimized substrate conditions . This implies that the calibrated respiration rate value lumps together maintenance (renewal of biomass synthesis) and growth (new biomass synthesis) respiration. We conducted simulations using the three different methods over a 5-day period. For the LBM-based method a time step of 0.43 seconds is used for diffusion and the synchronous transformation processes with the same time step . For the ball network model, we employed the calibrated implicit diffusion scheme from  and the asynchronous algorithm for transformation simulation with a 5-second time step. For our approach, we conducted tests using different time steps to understand their impact on simulation results (see Appendix S2 in S1 File for more details). Table S2T1 provides details on these time steps along with the elapsed time for each simulation. On a standard PC with an AMD Ryzen 7 PRO 6850H processor and 32.0 GB of RAM, the computation times are approximately as follows: LBM-based simulation: 3 weeks VGA simulation using an explicit scheme (test 3): 7.3 days VGA simulation using an implicit scheme (test 1): 5.2 days (refer to Table S2T1 in Appendix S2 in S1 File for more details) Pore Network Model (PNGM)-based method: 45 minutes The optimal intercorrelation was achieved in test 3 (see Appendix S2 in S1 File for more details), which is expected because it used the same time step as the LBM-based simulation, specifically 0.43 seconds for transformation. The decision to use a time step of 0.1 seconds for the explicit scheme was based on two factors: first, the accuracy of the explicit scheme, as it does not require approximation; second, 0.1 seconds represents the largest possible time step for the explicit scheme (at least for this dataset and experimental setup), otherwise we get negative values. In contrast, Test 4 demonstrates that an implicit scheme with the same time step as that used for the pore network model (5s) produces results that are significantly divergent from both LBM and PNGM simulations (Fig S2C in S1 File). The voxel graph-based approach presents a practical alternative to LBM for simulation of microbial decomposition of organic matter in soil, as it yields comparable results at lower computational costs. The results obtained indicate that VGA outperforms the other two methods by striking a balance between accuracy and computational time. This superiority can be attributed to three main factors: firstly, the representation of the pore space using a graph of connected voxels, employing only a 6-connectivity policy, results in a more accurate representation than, the pore network models and a lighter representation compared to the lattices used in LBM (D3Q19 in this case). Secondly, simulating mass transport between connected voxels using Fick’s law is computationally less expensive than simulating discrete collisions and particle propagation within a lattice grid. Finally, the implicit scheme of the graph diffusion equation, which allows the use of larger time steps than those used in LBM, is a key factor in reducing the number of iterations. The PNGM-based approaches demonstrate significant cost-effectiveness compared to voxel-based approaches such as the LBM-based approach or the discussed voxel graph-based approach. Even time steps of 30 seconds still yield good results (refer to  for more details). In the next section, a detailed theoretical approach will be provided for calculating the diffusional conductance coefficients for PNGM models using the voxel graph-based approach, in order to improve accuracy. 3.3. Improving PNGM simulations In the context of the ball network detailed in Section 2.4, we initialize the parameters denoted as Θ = θ i , j : i , j ∈ E P N G M by ∀ i , j ∈ E P N G M , θ i , j = S i , j d i , j Where S i , j represents the contact surface and d i , j signifies the distance between the two connected primitives i and j . Subsequently, a synthetic 3D pore space image is generated to represent the balls, followed by the construction of the voxel graph as outlined in section 2. Thirty scenarios of Distribution of Organic Matter (DOM) within the pore space voxels are generated as follows: for each scenario, a random mass M 0 is generated and distributed randomly among the voxels. The diffusion process is simulated using the explicit scheme of the voxel graph diffusion equation with a time step of 0.1 seconds and a diffusion coefficient D c = 100950 v o x e l 2 . j − 1 . The mass distribution in the voxels is recorded every 10 seconds during the simulation. Subsequently, all resulting voxel distributions are mapped to the corresponding ball distributions. These ball distributions are used to minimize the implicit objective L 2 (Eq. ). A total of 3600 data points is obtained, each consisting of a mass distribution in the ball network and the corresponding mass distribution after 10 seconds of diffusion simulated using the accurate VGA-based simulation. We run 1000 epochs of stochastic gradient descent for training. Each epoch involves the selection of four random data points from the obtained dataset, then calculating the normalized gradient and updating the learning parameters accordingly. Training starts with a learning rate of 0.1, which is then halved every 10 epochs. During training, we calculate the L 2 loss applied to four chosen data points to track minimization. The history of the L 2 values is plotted in . The loss decreases significantly indicating rapid learning and convergence within the first 10 epochs. Also, the training shows great stability and convergence to the solution, with a minimal L 2 value of 1.37766510 × 10 − 8 . After epoch 355, the learning process stabilizes, and the stochastic gradient descent algorithm converges to the set Θ = θ i , j : i , j ∈ E P N G M that minimize the application Θ ↦ L 2 Θ for the data provided. The obtained coefficients replace the coefficients α i , j : i , j ∈ E P N G M in one of the schemes (either Eq. or Eq. ) to simulate diffusion within the ball network, which serves as an approximation of the pore space. Since the number of spheres is considerably smaller than the number of voxels in the pore space, simulations using spheres require less computational effort than those using voxels. However, accurate diffusion simulations within the sphere network rely on precisely determining the diffusional conductance coefficients . and compare simulation results using two different sets of diffusional coefficients: the previous coefficients from  and the new coefficients derived using stochastic gradient descent based on data generated from the voxel approach. These results are then compared with the results of the accurate VGA. presents a comparison of simulated mass after 1.73 hours of diffusion, using the same scenario described in Section 3.1, with both VGA-based and PNGM-based simulations. The figure includes results from VGA using explicit scheme with a 0.1s time step a PNGM model calibrated by LBM simulation (referred to as the old PNGM) and from a PNGM model with diffusional conductance coefficients obtained after 1000 epochs of training (referred to as the new PNGM) using explicit scheme with a 0.1s time step. The new PNGM shows an improvement over the old PNGM, demonstrating enhanced accuracy in the simulation results. The difference between the VGA and PNGM curves near the origin is due to the limited number of spheres intersecting the first plane, causing the initial DOM mass to be more concentrated in the second plane. Specifically, in the VGA approach, for each voxel containing organic matter, the nearest sphere in the PNGM is identified, and the corresponding organic matter is placed accordingly. We conducted microbial decomposition simulations under the same scenario as in Section 3.2, employing the old PNGM-based simulation calibrated using LBM in previous work , and the improved PNGM-based simulation using the newly obtained diffusional conductance coefficients . These simulations were done using the implicit scheme with a 10s-time step and the asynchronous transformation using a 30s-time step. The results of these simulations are compared with those of the VGA-based simulation using the explicit scheme with a 0.1s-time step and the asynchronous scheme of transformation with a 0.43s-time step. compares simulated microbial decomposition, using both VGA and PNGM. The figure includes results from PNGM calibrated using LBM simulations as described before (labeled “old PNGM”), PNGM with diffusional conductance coefficients optimized through training (labeled “new PNGM”), and VGA. While the “new PNGM” model maintains the structural characteristics of the “old PNGM” model, its curves are consistently closer to the VGA-based simulation results across all components of microbial decomposition. This indicates that the optimization of diffusional conductance coefficients improved the model’s performance. The following table summarizes the pros and cons of all the simulations approaches discussed in this paper. 

To compare diffusion simulation using the graph diffusion equation with simulations using the Lattice Boltzmann method and the pore network model, we draw upon a prior experience detailed in , which was conducted to compare diffusion using a pore network model (ball network model) and LBM, and we outline our approach as follows: Uniformly, we distributed M 0 = 592.7593 m g of DOM on the first two z layers ( z = 1 and z = 2 ) of the 512 × 512 × 512 image. In our approach and the LBM, we selected all the voxels of the two layers ( z = 0 and z = 1 ) corresponding to the pore space and distributed M 0 uniformly among the selected voxels. For the ball network model, we selected the balls intersecting the two layers and distributed the mass M 0 homogeneously across them (i.e., applying the same concentration to all the balls). The experiment aims to compare mass profiles in each layer of the image after some time of 1.76 h o u r s . We ran the simulations using the three methods and obtained the final distributions, from which we calculated the mass in each layer of the image. Subsequently, we plot these results together for comparison. The basic principle of pore network geometrical modelling (PNGM) is to globally transport mass from one ball to the connected balls following Fick’s laws. Therefore, we need to calibrate the diffusional conductance. We employed a general coefficient α = 0.6 , which we multiplied by the contact surface area between all pair of geometrically connected balls, calibrated from LBM . Details on how the coefficient α is determined are provided in Section 2.4. illustrates a global comparison of the overall accuracy of each method. PNGM was simulated using the implicit scheme with a time step of 15 seconds, while GDE was solved using the implicit scheme with a time step of 30 seconds. In contrast, the time step for LBM is 0.43 seconds, calculated based on the 3D image resolution. The results obtained with the implicit scheme of the GDE are very close to those of LBM and much more accurate than those of the implicit scheme of PNGM. The use of voxels to discretize the space is similar to LBM; however, our approach models diffusion as mass transport according to Fick’s law, whereas LBM uses particle-based collision and propagation processes. Although both methods are conceptually similar, our approach involves less calculations. In contrast, PNGM uses geometrical primitives that are larger than a single voxel and vary in size. This approach makes modeling mass transport in the pore space using Fick’s law less precise. The explicit scheme for solving the graph diffusion equation requires the use of small-time steps to maintain stability and accuracy in the numerical solution. This necessity arises because each time step needs to be small enough to ensure that the diffusion process is accurately captured and does not lead to instability or inaccuracies in the results. As a result, the explicit scheme involves frequent updates and calculations, which can become computationally expensive, especially for large-scale problems or long simulation times. In contrast, the implicit scheme offers a more computationally efficient alternative by allowing the use of extended time steps. This approach improves computational efficiency because it reduces the number of required updates and calculations, as larger time steps cover a greater extent of the diffusion process in each iteration. Although implicit methods are typically more complex to implement, they often provide better stability and accuracy over longer time periods. illustrates the comparison between these two schemes. In the explicit scheme, a time step of 0.1 seconds is used, resulting in frequent updates and higher computational costs. On the other hand, the implicit scheme employs a relatively extended time step of 30 seconds, demonstrating a reduction in computational effort while still maintaining accuracy in the simulation results. 

In this section, we validate the efficiency of our method for simulating microbial decomposition of organic matter by comparing its results with those obtained using the LBM-based approach and the PNGM-based method. We conducted numerical simulations of microbial decomposition using the LBM-based approach, the Voxel Graph-based Approach (VGA), and the ball network-based approach on the same dataset previously discussed. The simulation procedure is as follows: Initially, we evenly distributed M 2 0 = 289.5 μ g of DOM and introduced 1000 bacterial spots (clusters) into the pore space at random locations. This patchiness of bacterial distribution, as a result of cell growth mechanisms and environmental constraints, has been investigated in . The clusters collectively represent 5.2 × 10 7 bacteria cells, which correspond to M 1 0 = 5.2 × 10 7 × 5.41 × 10 − 8 μ g = 2.8132 μ g of carbon. For the PNGM-based approach, we distributed the mass M 2 0 homogeneously, (i.e., the same concentration in each ball), and we put bacterial mass into the balls corresponding to the random location of the spots. The total mass in the pore space initially is as follows: M 1 0 , M 2 0 , M 3 0 , M 4 0 , M 5 0 , where M 3 0 = M 4 0 = M 5 0 = 0 , corresponding to 99.0376 % of dissolved organic matter and 0.962 % of living microorganisms. The biological parameters employed in this study were adopted from  for Arthrobacter sp. 9R that was calibrated from real experimental data. These parameters are as follows: ρ = 0.2 d a y − 1 is the relative respiration rate, μ = 0.5 d a y − 1 is the relative mortality rate, β = 0.55 is the proportion of microbial mass (MB) returning to DOM (the remaining portion is for MB returning to SOM), v S O M = 0.01 d a y − 1 and v F O M = 0.3 d a y − 1 are the relative decomposition rates for SOM and FOM respectively, v D O M = 9.6 d a y − 1 and K D O M = 0.001 g C . g − 1 are the maximum relative growth rate of MB, and the half saturation constant of DOM respectively. D D O M = 100950 v o x e l 2 . d a y − 1 is the diffusion coefficient of DOM in water. The value used corresponds to the diffusion coefficient of an organic molecule in water . Parameterization of our model has been performed on experimental measurements of growth and respiration curves of Arthrobacter sp cultivated in optimized substrate conditions . This implies that the calibrated respiration rate value lumps together maintenance (renewal of biomass synthesis) and growth (new biomass synthesis) respiration. We conducted simulations using the three different methods over a 5-day period. For the LBM-based method a time step of 0.43 seconds is used for diffusion and the synchronous transformation processes with the same time step . For the ball network model, we employed the calibrated implicit diffusion scheme from  and the asynchronous algorithm for transformation simulation with a 5-second time step. For our approach, we conducted tests using different time steps to understand their impact on simulation results (see Appendix S2 in S1 File for more details). Table S2T1 provides details on these time steps along with the elapsed time for each simulation. On a standard PC with an AMD Ryzen 7 PRO 6850H processor and 32.0 GB of RAM, the computation times are approximately as follows: LBM-based simulation: 3 weeks VGA simulation using an explicit scheme (test 3): 7.3 days VGA simulation using an implicit scheme (test 1): 5.2 days (refer to Table S2T1 in Appendix S2 in S1 File for more details) Pore Network Model (PNGM)-based method: 45 minutes The optimal intercorrelation was achieved in test 3 (see Appendix S2 in S1 File for more details), which is expected because it used the same time step as the LBM-based simulation, specifically 0.43 seconds for transformation. The decision to use a time step of 0.1 seconds for the explicit scheme was based on two factors: first, the accuracy of the explicit scheme, as it does not require approximation; second, 0.1 seconds represents the largest possible time step for the explicit scheme (at least for this dataset and experimental setup), otherwise we get negative values. In contrast, Test 4 demonstrates that an implicit scheme with the same time step as that used for the pore network model (5s) produces results that are significantly divergent from both LBM and PNGM simulations (Fig S2C in S1 File). The voxel graph-based approach presents a practical alternative to LBM for simulation of microbial decomposition of organic matter in soil, as it yields comparable results at lower computational costs. The results obtained indicate that VGA outperforms the other two methods by striking a balance between accuracy and computational time. This superiority can be attributed to three main factors: firstly, the representation of the pore space using a graph of connected voxels, employing only a 6-connectivity policy, results in a more accurate representation than, the pore network models and a lighter representation compared to the lattices used in LBM (D3Q19 in this case). Secondly, simulating mass transport between connected voxels using Fick’s law is computationally less expensive than simulating discrete collisions and particle propagation within a lattice grid. Finally, the implicit scheme of the graph diffusion equation, which allows the use of larger time steps than those used in LBM, is a key factor in reducing the number of iterations. The PNGM-based approaches demonstrate significant cost-effectiveness compared to voxel-based approaches such as the LBM-based approach or the discussed voxel graph-based approach. Even time steps of 30 seconds still yield good results (refer to  for more details). In the next section, a detailed theoretical approach will be provided for calculating the diffusional conductance coefficients for PNGM models using the voxel graph-based approach, in order to improve accuracy. 

In the context of the ball network detailed in Section 2.4, we initialize the parameters denoted as Θ = θ i , j : i , j ∈ E P N G M by ∀ i , j ∈ E P N G M , θ i , j = S i , j d i , j Where S i , j represents the contact surface and d i , j signifies the distance between the two connected primitives i and j . Subsequently, a synthetic 3D pore space image is generated to represent the balls, followed by the construction of the voxel graph as outlined in section 2. Thirty scenarios of Distribution of Organic Matter (DOM) within the pore space voxels are generated as follows: for each scenario, a random mass M 0 is generated and distributed randomly among the voxels. The diffusion process is simulated using the explicit scheme of the voxel graph diffusion equation with a time step of 0.1 seconds and a diffusion coefficient D c = 100950 v o x e l 2 . j − 1 . The mass distribution in the voxels is recorded every 10 seconds during the simulation. Subsequently, all resulting voxel distributions are mapped to the corresponding ball distributions. These ball distributions are used to minimize the implicit objective L 2 (Eq. ). A total of 3600 data points is obtained, each consisting of a mass distribution in the ball network and the corresponding mass distribution after 10 seconds of diffusion simulated using the accurate VGA-based simulation. We run 1000 epochs of stochastic gradient descent for training. Each epoch involves the selection of four random data points from the obtained dataset, then calculating the normalized gradient and updating the learning parameters accordingly. Training starts with a learning rate of 0.1, which is then halved every 10 epochs. During training, we calculate the L 2 loss applied to four chosen data points to track minimization. The history of the L 2 values is plotted in . The loss decreases significantly indicating rapid learning and convergence within the first 10 epochs. Also, the training shows great stability and convergence to the solution, with a minimal L 2 value of 1.37766510 × 10 − 8 . After epoch 355, the learning process stabilizes, and the stochastic gradient descent algorithm converges to the set Θ = θ i , j : i , j ∈ E P N G M that minimize the application Θ ↦ L 2 Θ for the data provided. The obtained coefficients replace the coefficients α i , j : i , j ∈ E P N G M in one of the schemes (either Eq. or Eq. ) to simulate diffusion within the ball network, which serves as an approximation of the pore space. Since the number of spheres is considerably smaller than the number of voxels in the pore space, simulations using spheres require less computational effort than those using voxels. However, accurate diffusion simulations within the sphere network rely on precisely determining the diffusional conductance coefficients . and compare simulation results using two different sets of diffusional coefficients: the previous coefficients from  and the new coefficients derived using stochastic gradient descent based on data generated from the voxel approach. These results are then compared with the results of the accurate VGA. presents a comparison of simulated mass after 1.73 hours of diffusion, using the same scenario described in Section 3.1, with both VGA-based and PNGM-based simulations. The figure includes results from VGA using explicit scheme with a 0.1s time step a PNGM model calibrated by LBM simulation (referred to as the old PNGM) and from a PNGM model with diffusional conductance coefficients obtained after 1000 epochs of training (referred to as the new PNGM) using explicit scheme with a 0.1s time step. The new PNGM shows an improvement over the old PNGM, demonstrating enhanced accuracy in the simulation results. The difference between the VGA and PNGM curves near the origin is due to the limited number of spheres intersecting the first plane, causing the initial DOM mass to be more concentrated in the second plane. Specifically, in the VGA approach, for each voxel containing organic matter, the nearest sphere in the PNGM is identified, and the corresponding organic matter is placed accordingly. We conducted microbial decomposition simulations under the same scenario as in Section 3.2, employing the old PNGM-based simulation calibrated using LBM in previous work , and the improved PNGM-based simulation using the newly obtained diffusional conductance coefficients . These simulations were done using the implicit scheme with a 10s-time step and the asynchronous transformation using a 30s-time step. The results of these simulations are compared with those of the VGA-based simulation using the explicit scheme with a 0.1s-time step and the asynchronous scheme of transformation with a 0.43s-time step. compares simulated microbial decomposition, using both VGA and PNGM. The figure includes results from PNGM calibrated using LBM simulations as described before (labeled “old PNGM”), PNGM with diffusional conductance coefficients optimized through training (labeled “new PNGM”), and VGA. While the “new PNGM” model maintains the structural characteristics of the “old PNGM” model, its curves are consistently closer to the VGA-based simulation results across all components of microbial decomposition. This indicates that the optimization of diffusional conductance coefficients improved the model’s performance. The following table summarizes the pros and cons of all the simulations approaches discussed in this paper. 

This study introduces the Voxel Graph-based Approach (VGA) to simulate microbial activity using the set of voxels of pore spaces obtained using 3D computed tomography images of soil. Traditional methods such as the Lattice Boltzmann Method (LBM) or PDE solvers are computationally demanding when dealing with complex 3D structures or large image data. Pore Network Geometrical Models (PNGM), while more efficient, often require calibration and suffer from reduced accuracy due to oversimplified geometrical representations. VGA addresses these issues by considering an exact representation of the pore space using a graph of connected voxels directly from 3D micro-CT images, allowing accurate modeling of transport processes based on Fick’s law. When coupling the transport processes with biodegradation model, our approach simulates microbial decomposition of organic matter in saturated porous media accurately and faster than traditional methods. Our comparative analysis demonstrates that VGA achieves comparable accuracy to LBM while reducing computation time by a quarter. Although VGA is slower than PNGM, it provides superior accuracy and requires no calibration. PNGM-based simulations can yield results close to those of VGA or LBM when a robust calibration method is applied. Notably, we showed that VGA improves PNGM-based simulations by refining diffusional conductance coefficients through stochastic gradient descent, which can be used to improves the precision of pore network geometrical modeling in large-scale applications. The VGA approach offers several key advantages: Refining Diffusional Conductance Coefficients: VGA generates high-accuracy data with minimal computational cost, which can be used to approximate diffusional conductance coefficients in PNGM-based simulations. Incorporating Pore Space Deformation: Unlike PNGM, which relies on rigid geometric primitives, VGA represents pore space as a graph of connected voxels. This representation allows for dynamic changes, such as voxel movement or removal, enabling the modeling of pore space deformation. Upscaling: By integrating VGA with graph neural networks, we can train models to predict long-term microbial dynamics from 3D images of larger samples, addressing upscaling challenges effectively. Generalizing: VGA can simulate a wide range of diffusion and transformation processes in two-phase porous media or fractured environments. However, VGA also has limitations: Spatial Discretization: The accuracy of the model is influenced by the spatial discretization, determined by the voxel resolution. In our case, the voxel size of 24 µm may limit the resolution of small-scale features. Two-Phase System: The current implementation is restricted to a two-phase system, which may not fully capture the complexity of multi-phase systems found in real-world scenarios. However, it is straightforward to extend the framework by incorporating additional phases using a heterogeneous graph approach, which we plan to explore in the future. Advective Transport: The VGA framework does not account for advective transport processes, which can be significant in many environmental and biological systems. Despite the limitations, VGA provides a valuable tool for simulating microbial decomposition of organic matter in complex 3D soil structures. Compared to recent methods by Pot et al. (2022), Monga et al. (2022), and Zech et al. (2022) that is restricted to 2D applications, the voxel graph-based approach discussed in this paper is more efficient for simulating microbial decomposition of organic matter in 3D porous media obtained from computed tomography images of soil. 

S1 File Appendices S1, S2 and S3. (DOCX) 

