Modeling Microbial Decomposition in Real 3 D Soil Structures Using Partial Differential Equations

Partial Differential Equations (PDEs) have been already widely used to simulate various complex phenomena in porous media. This paper is one of the first attempts to apply PDEs for simulating in real 3D structures. We apply this scheme to the specific case study of the microbial decomposition of organic matter in soil pore space. We got a 3D geometrical representation of the pore space relating to a network of volume primitives. A mesh of the pore space is then created by using the network. PDEs system is solved by free finite elements solver Freefem3d in the particular mesh. We validate our PDEs model to experimental data with 3D Computed Tomography (CT) images of soil samples. Regarding the current state of art on soil organic matter decay models, our approach allows taking into account precise 3D spatialization of the decomposition process by a pore space geometry description.


Introduction
It is a fact that soil structure is a heterogeneous media and organic matter decomposition in soil is one of the most complex ecological processes.The nature of organic matter, the dynamic of microorganisms, the environmental conditions etc. are the major factors which influence the decomposition of organic matter.Despite that microbial dynamics plays an important role in soil organic matter decomposition, the majority of organic matter models make the assumption that carbon limitation is controlled only by the intrinsic degradability of organic matter [1].There are few models that take into account the rule of exoenzymes produced by microorganisms to convert complex substrates into available compounds [2].Moreover, most of organic matter models do not consider the physical heterogeneity that controls partly the availability of organic matter to microorganisms [3].
In the recent contributions [4,5], the authors tried to simulate the diffusion of soluble carbon substrate in soil in 2D or 1D.Some attempts were also made to simulate enzyme diffusion in artificial structured environments [6,7].But none of these models consider explicitly the real 3D soil structure.It is obvious that soil organic matter modeling needs to be in 3D soil structure because 1D or 2D representations of soil structure [8,9] can only supply partial information on the nature of overlapping [10].Yet the improved performance of 3D X-ray computed tomography sensors now makes it possible to obtain very high resolution images of soil sample volume.Indeed, such non-destructive imaging technologies (Xray tomography or CAT, laser scanning confocal microscopy or LSCM) can generate digital gray scale images that visualize 3D soil structure at aggregate [11] and also smaller intra-aggregate scale (16 μm in [12]; 3 μm in [13]).
Recently, the new model of Monga et al. in [14] simulated the decomposition of organic matter at a microscale using computer tomography images of real 3D soil but it does not take into account the diffusion process of dissolved organic matters.The main reason is that it costs a lot of time for computing diffusion process in 3D soil structure even with a very strong computer.
In this study we present an organic matter model that focuses on both transformation and diffusion processes in a 3D soil structure.We use Partial Differential Equations (PDEs) in order to simulate biological activity.In fact, models available in literature that simulates the biological activity in porous media use different techniques like cellular automata to simulate biomass growth in biofilms [15], Partial Differential Equations (PDEs) to simulate biomass recycling in fungal colonies [16] and multiagent system to simulate the impact of earthworms on soil structure [17].In this work, we focus on the PDE method because 1) PDEs have been widely used to simulate various complex phenomena, particularly including transformation and diffusion processes in porous media but few studies have used this approach to biological soil systems; 2) the recently increasing performances of computer PDE solvers make it possible to simulate more and more complex systems.
This paper dealt with the application of PDEs to simulate organic matter degradation in real micro 3D structures and also took the diffusion into account.We tested our method using high resolution 3D CT of real soil images.Our work faces the problem of solving a nonlinear PDE system (reaction-diffusion) in a non-regular 3D mesh.
The paper is organized as follows.Section 2 presents the biological model.Section 3 presents the mathematical model of soil organic matter decay using PDE system.The PDE system is formed by reaction diffusion equations.The approximation of the model weak solutions is found by using a finite element method and a Newton algorithm.In Section 4, we implement the resolution algorithm by using Freefem3d software [18,19].In Section 5, we present a validation of our PDE model to a experimental result.Section 6 presents a conclusion and some perspectives.Detailed Freefem3d code is given in Appendix.

Biological Model
Our aim is to simulate microbial decomposition in a nonregular 3D geometric volume shape defined by pore space.In this case pore space is described by a set of volume primitives [14,[20][21][22].Therefore, we take as input the following data:  Geometrical representation of pore space using a network (attributed relational valuated graph) volume primitives [14,[20][21][22]. Parameters describing the initial spatial distribution of elements involved in biological activity: microbial biomass (MB), fructose (F), solid microbial wastes (SW), decomposed soluble microbial metabolites (DSM), un-decomposed soluble microbial metabolites (UDSM) and inorganic organic matter (CO 2 ).Thus, we assume that the microbial decomposition process involves six biological elements noted MB, F, SW, DSM, UDSM and CO 2 .We denote F the amount of fructose which is added into the environment.And we assume that it diffuses through water path to be consumed by microbial biomass (MB) for their maintenance and growth.It is supposed that MB does not move, so we assume that their diffusion coefficient is very small comparing with others.We further assume that the dead MB is transformed into three kinds of organic matters: solid microbial wastes (SWs) which cannot be consumed by MB and cannot diffuse, un-decomposed soluble microbial metabolites (UDSMs) which can diffuse but cannot be consumed by MB, decomposed soluble microbial metabolites (DSMs) which can be consumed again by MB and can diffuse.MB breathes to produce inorganic carbon (CO 2 ) [2]. Figure 1 summarizes the biological process inside one spatial unit.The output of the simulation system for each step time is the precise distribution of biological activity parameters, i.e.: MB density, F density, SW density, DSM density, UDSM density and CO 2 density.
Therefore, we provide a kind of animated film showing spatially the evolution of biological characteristics.From these information, we can easily compute the classical global evolution curves: CO 2 content, UDSM content, MB content...    be a point of the pore space.We

Presentation of the Model Using PDE
and   , c x t the densities of microbial (MB), of fructose (F), of solid microbial wastes (SW), of un-decomposed soluble microbial metabolites (UDSM), of decomposed soluble microbial metabolites (DSM) and of inorganic organize matter (CO 2 ) at position and at time , respectively.
x   t In the next subsections, we presented replicator equations for the variations of quantities of the six biological elements.

Microbial Biomass (MB)
Let be an elementary volume in .The variation of the quantity of MB in is due to: microbial diffusion, microbial growth, microbial mortality and microbial breathing.Thus during the breathing, microbial decomposers lose a part of their masses.The following equation summarizes the process: We assume that the MB growth depends on the quantities of fructose (F) and of decomposed soluble microbial metabolites (DSMs).We also assume that microbial consumes F and DSMs according to the Monod equation, proposed in J.D. Murray's book [23], as follows: where the variables are set as follows.b is the concentration of microbial decomposers, k is the maximal growth rate, b K is the half-saturation constant, f is the concentration of fructose, n is the concentration of decomposed soluble microbial metabolites.
The above equation describes microbial decomposers growth rate.Thus, the variation of b is expressed by the following equation: where is the diffusion operator, b represents the microbial decomposers diffusion coefficient, the mortality rate, is set to the breathing rate.r

Fructose (F)
F density variation comes from: F diffusion and con-following equation summarizes F density variation process: F variation is ruled by the following equation: , where is the diffusion operator, where is the rate of dead MB trans-

Metabolites (UDSM)
M quantity variation is ca of a part of MB mortality and the diffusion process as follows:  UDSM quantity variation is expressed by the following equation: sumption of F by microbial decomposers.Thus, the where is the diffusion operator,  is the rate of de ran in ad MB t sformed DSM ation comes from: DSM diffusion, a to UDSM.

Decomposed Soluble Microbial Metabolites (DSM)
density vari part of dead MB and consumption of DSM by microbial decomposers.Thus, the following equation summarizes DSM density variation process: DSM variation is ruled by the following equation: , where Δ is the diffusion operator, represents DSM n D diffusion coefficient and  is th ate of dead MB transformed into DSM.

Inorganic Carbon (CO e r
2 ) e to its diffusion Inorganic carbon (CO 2 ) variation is du and to its production by microbial decomposers (MB) during breathing.
The CO 2 dynamics equation reads as follows: where is the diffusion operator, is set to the ivities act only insde

Boundary Conditions
We assume that biological act  he and there is not any outside forces that effect on t dynamics inside  .It means that flow is null on  for all variables.erefore, on the border of Th  no  , we use the Neumann boundary conditions.

The Complete P uations describing
In this section we use the above eq variable variations to set a global PDE system modeling biological dynamics.Let 0 T  be a fixed time and let's define Therefore, the whole system of partial differential equations governing the biological model becomes in T  : , , .
We use Neumann hom and the following initial conditions in (14) ogeneous boundary conditions he functions of or CO 2 .The in position x , i.e. amount of the bi ts at the beginning are known at any position in ological elemen  .Therefore, the above PDE system describes precisely microbial decomposition o organic matter in soil.In the fo f g it into a g way: The diffusion coefficients matrix llowing section, we show how to solve this PDE system which can practically simulate soil biological activity.


The reaction terms of equations are represented functions defined as follows: by , 1,2, ,6     Let's define the vector function F such that The vector form of the system is

lation
Let's introduce the following Sobolev space in , Assuming that data are sufficiently regular, the variational formulation consists in finding a function (26)  

Building a Mesh h  of Domain 
A mesh of pore space    can be built from a georepresentatio metrical n pore space [14,[20][21][22].The following algorithm can then be use to build a mes  Computing the geometrical primitives netwo presenting pore space jac resenting pore ere  of h: rk re-We follow the approach described in [14,[20][21][22].In detail, we compute an ad ency graph rep space using 3D computed tomography image of a soil sample.Geometrical primitives can be chose, for instance, the balls as in [12].Therefore, we obtain an adjacency graph of balls describing pore space wh each node is attached to a ball representing a pore within common sense [14,[20][21][22].
Building a mesh by using all centers of the balls From all centers of the balls, we build a tetrahedra mesh for  .In order to do that, we propose using the Delaunay triangulation program developed at INRIA by the GAMMA project [24,25].
Building the mesh following this way is exactly the same idea as in [14,22] where the transformation process acts inside each ball and the diffusion process acts trical m er el biological activities. between the connected balls.Instead of the geome odelling approach using updating graph in [14,22], here we consid in this work the mathamatical modelling approach.We are trying to use the benefits of the approach in order to better mod Pruning the mesh We do pruning the mesh in order to neglect edges between two centers of balls that do not connect.Now we get a mesh h  of domain  .After building a mesh h  of domain  , we solve the variational formulation (27) in the following discrete space: .
n the finite dimensional space .The problem consists in solving the following

Numerical Scheme
Th of e numerical resolution he problem is ivided into three steps:  Step 1: We discretize the problem via finite element method i , , and We consider the following approximation: The numerical scheme becomes:

Resolution Algorithm
Fr en in order to itialize to .Afterwa s, we repeat th om the initial conditions, we know U .Th find solution n U knowing of the loop for linearization.

8)
a When one of the stopping criteria is satisfied, * U represents the se hed solution hus, the resoluti algo m can be described as follows: 1) max N and a  re fix 2) If the cond r is satisfied, the lo , we give the Freefem3d code we implement for solving the system.Figure 2 shows slices (1650 μm × 1650 μm) of the 3D volume image.We extract from the (1650 μm × 1650 μm × 1745 μm) image a (400 μm × 400 μm × 400 μm) image for memory requirements and computing time reasons.Figure 3 shows successive slices (400 μm × 400 μm) of the (400 μm × 400 μm × 400 μm) 3D image extracted.

Experimental Results on Real Data
The second stage provides a set of voxels corresponding to pore space.The threshold is computed using the expected porosity of the sample.Figure 4 shows a cross section representing pore space (white color), the porosity is 35%.Using the geometric model for the 3D representation of the pore space described in [14,[20][21][22], we calculate the network of balls.
Afterwards the balls are drained for the following values of water pressure ψ = −0.01MPa and ψ = −0.001Mpa corresponding respectively to situations that are relatively wet (high water pressure) and dry (low water pressure), respectively.For water pressure ψ = −0.01MPa we get 649776 balls and for water pressure ψ = −0.001Mpa we get 569175 balls.Figure 5 shows perspective views of the ball based pore space representation.Figure 6 shows the tetrahedral mesh building from all the centers of the balls.
We built the model for two different values of water pressures (high water pressure, wet and low water pressure, dry).Initial microbial biomass (40 g) is put on the mesh such that probability of microbial being in a node is proportional to the radius of the bal which corresponds to the node.At the beginning, 300 g fructose is uniformly distribute in the mesh.Diffusion coefficients of fructose, CO 2 and DSM are equal to  which is very small in comparison with diffusion coefficient of the others since we assume that microbial does not move.We compare experimental carbon flow curves as well as dissolved organic matter (DOC) amount with the ones provided by our PDE model.We test our PDE model for five different types of microbial species: Arthrobacter sp 3R, 7R and 9R, Rhodococcus sp 6L and 5L.Parameter values are taken from the investigation of Coucheney in [27] r given in Table 1. Figure 7          experimental data in the wet condition (left), while grey solid with circle marks and tetra-angle marks are about CO 2 and DOC from experimental data in dry condition (right).Similar comparison can be seen in Figure 9 for Rhodococcus sp 6L and 5L.It is noticed that we assumed that the sum of decomposable and undecomposable bacteria metabolites (DMM + UDMM) is simulated at the end of the experiment in order to measure the bound cases except the case Arthrobacter 3R and 9R with high water pressure.It can be seen that CO 2 fits very well to CO 2 -data for all cases.It is sligtly different between the two with Arthrobacter 7R in low water pressure.

Conclusion and Perspectives
In this study we have proposed a new model of organic matter decomposition in the 3D soil pores.The novelties of our approach are to consider the real 3D microstructures of soil, to take into account diffusion and transformation processes.To our knowledge, this is the first attempt to describe diffusion process in 3D soil structure thanks to the strength of PDEs.We applied our PDE model to real CT images of soil in which we add realistic amount of organic matter and microorganisms.We validated our model to experimental results.In particular, we compare experimental carbon flow curves as well as DOC amount with the ones provided by our PDE model at two different water pressures.The simulation results show that our PDE model fits reasonably well to the experimental results.In this current contribution, we just consider the dynamics of single microbial species.In future studies, simulations with real experimental data will be carried out with our model in order to analyze microbial competition of degradation nder different

Appendix. Algorithm Code Using Freefem3D
The freefem3d code corresponding to the algorithm is the following: Step 1.Getting the mesh from the source file.meshM = read(medit, "tetranuage.mesh") Step 2. Defining the data.Functions , 1,2, ,6   of type femfunction defined on the mesh i f i M are used to represent reaction terms of the system.To solve the nonlinearity of the reaction functions, we use the same type to define approximate functions eb, ef, em 1 , em 2 , en, and ec.The same function type is used to build terms correction db, df, dm 1 , dm 2 , dn, and dc.The pa of rameters ared M de to imple-such as diffusion coefficients, growth and death rate of MB are declared using double type.
Step 3. Defining and initializing the functions representing variable densities.Initial functions are decl with the femfunction type defined on M. The ofstream function permits to create a file to save solutions.The quantities of the variables are obtained by integrating their density in the pore space represented by the mesh Thus, we use the following Freefem3d co .ment the numerical scheme: // opening of a file named "outputpathname" to save results of stream out put = of stream("outputpathname"); // saving of the time and initial quantities  ; // loop for linearizing the system do { //we redefine the reaction terms as functions of approximate terms.// solving the system to find correction terms.// Db, Df, Dn, Dc are diffusion coefficients.
comes from the tr a part of MB mortality.1 variation of m  a part of mortality of .b Thus SW quantity variation is expressed b part of mortality of .
a maximal number of iterations  we define a very small positive real and we stop each ration.Two stopping criteria are de called max N ; We have validated our PDE mo Th ag RA eriment 50 km west of Paris) in France.The soil columns were scanned by means of a high resolution m and) operating at 90 KeV and a current of 112 mA. on real data of sand.e soil is sampled in the surface layer (0 -20 cm) of an ricultural field at the IN exp al site of Feu-cherolles (FEU) ( icro-CT machine (SIMCT Equipment: SIMBIOS, University of Abertay Dundee, Scotl Using X-ray tomography we compute a 3D Co Tomography (CT) image of a sand sample.The sizes of the original image are: 1650 μm × 1650 μm × 1745 μm.The voxel is cubic of resolution: 5 μm × 5 μm × 5 μm.

Figure 5 .
Figure 5. Perspective views of the ball based pore space representation.We display only the balls whose radius is at least 10.

Figure 6
Figure 6.A tetrahedral mesh building from the centers of 1000 balls.Table 1. Parameter values are used in the simulation where   k

Figure 7 .
Figure 7. Simulation results of four biological elements of Rhodococcus sp 6L are presented.The above figure is about simulation with high water pressure and the below figure is about simulation with low water pressure.In the above figure: grey solid with square marks is about experimental CO 2 , grey solid with tri-angle marks is about experimental DOC.In the below figure: grey solid with circle marks is about experimental CO 2 , grey solid with tetra-angle marks is about experimental DOC.Both: black solid is about CO 2 from the model, black dash is about DOC from the model, grey dark solid is about microbial biomass, grey dark dots is about fructose.

Figure 8 .
Figure 8.Comparison between model and experimental data for Arthrobacter sp 3R, 7R and 9R.On the left hand side is about high water pressure, on the right hand side is about low water pressure.Solid line is about CO 2 from the PDE mode, dash line is about DOC from the model, square is about experimental CO 2 , circle is about experimental DOC.

Figure 9 .
Figure 9.Comparison between model and experimental data forfor Rhodococcus sp 6L and 5L.On the left hand side is about high water pressure, on the right hand side is about low water pressure.Solid line is about CO 2 from the PDE mode, dash line is about DOC from the model, square is about experimental CO 2 , circle is about experimental DOC.
It was completely decomposed around the second day.The microbial biomass is therefore grows fast in the first two days.It is then decreases the days after.Since growth of CO 2 is propotional to the microbial biom to what we observe from the data of CO 2 fructose) are presented for Rhodococcus sp 6L.One can see in the figure that fructose was rapidly decomposed by microbial for both cases (high and low water pressures).