Simulating Biological Dynamics Using Partial Differential Equations: Application to Decomposition of Organic Matter in 3D Soil Structure

Abstract

The majority of carbon on earth is in the form of soil organic matter. And its degradation by microorganisms leads to the remineralization of carbon as carbon dioxide. The microbial activity causes a reduction of soil carbon and increases atmospheric carbon. However, most models of organic matter do not explicitly take into account this reality. We try to answer these questions by developing and validating a model describing the action of microorganisms on degradation of organic matter. We use simulation domain as the pores in the soil modeled by a network of balls. The model is solved numerically in the balls by the finite element method with the solver of partial differential equations (PDEs) Freefem3d. We compare the numerical results with experimental data on the mineralization of soil carbon.

Appendices

Appendix A: Numerical Scheme

1.1 A.1 Variational Formulation

We simplify the system writing by transforming it into a vector form by defining:

$$u\equiv\left( u_{1},u_{2},u_{3},u_{4},u_{5},u_{6}\right)^{t} \equiv\left( b,n,m_{1},m_{2},e,c\right)^{t} $$

and

$$u_{0}\equiv\left( b_{0},n_{0},m_{10},m_{20},e_{0},c_{0}\right)^{T}. $$

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

$$\underline{D}=\left( \begin{array}{cccccccc} D_{b} &0 &0 &0 &0 &0\\ 0 &D_{n} &0 &0 &0 &0\\ 0 &0 &0 &0 &0 &0\\ 0 &0 &0 &0 &0 &0\\ 0 &0 &0 &0 &D_{e} &0\\ 0 &0 &0 &0 &0 &D_{c}\\ \end{array} \right). $$

The reaction terms of equations are represented by functions F _i , i = 1,2,…,6 defined as follows:

$$\begin{array}{@{}rcl@{}} F_{1}(u)&=& \frac{ku_{1}}{K_{b}+u_{1}}u_{2} - (\mu+r+\nu)u_{1},\\ F_{2}(u)&=&\frac{ku_{5}}{K_{m}+u_{5}}\left( c_{1}m_{1}+c_{2}m_{2}\right)- \frac{ku_{1}}{K_{b}+u_{1}}u_{2}+\alpha_{1}(\zeta)u_{5}+\alpha_{2}(\mu)u_{1},\\ F_{3}(u)&=& -\frac{c_{1} u_{5}}{K_{m}+u_{5}}u_{3}+(1-\alpha_{1}(\zeta))u_{5}+(1-\alpha_{2}(\mu))u_{1},\\ F_{4}(u)&=& -\frac{c_{2} u_{5}}{K_{m}+u_{5}}u_{4},\\ F_{5}(u)&=& \nu u_{1}-\zeta u_{5},\\ F_{6}(u)&=& r u_{1} \end{array} $$

and

$$F(u) = \big(F_{1}(u),F_{2}(u),F_{3}(u),F_{4}(u),F_{5}(u),F_{6}(u)\big)^{T}. $$

The vector form of the system is

$$\left\{ \begin{array}{ll} \partial_{t}u =\text{div}\left( \underline{D}\nabla u\right) + F(u)&\text{ in}\,\,{\Omega}_{T},\\ \frac{\partial u}{\partial n}=0&\text{ on}\,\,\partial{\Omega}\times]0,T[,\\ u(t=0) =u_{0}&\text{ in}\,\,{\Omega}. \end{array} \right. $$

Assuming that data are sufficiently regular, the variational formulation consists in finding a function u(t)∈(H ¹(Ω))⁶ such that:

$${\int}_{\Omega}\frac{\partial u}{\partial t}vdx +{\int}_{\Omega}\underline{D}\nabla u\nabla vdx={\int}_{\Omega}F(u)vdx\qquad\forall v\in\left( H^{1}({\Omega})\right)^{6}. $$

After building a mesh Ω_h of domain Ω, we solve variational formulation in the following discrete space:

$$V_{\textit{h}}=\left\{v\in \left( C(\overline{\Omega})\right)^{6}: \forall K\in{\Omega}_{\textit{h}}\quad\left( v_{|K}\in P_{1}\right)^{6}\right\}. $$

A.2 Numerical Scheme

We discretize the problem via finite element method in the finite dimensional space V _h . The problem consists in solving the following system

$$\frac{\partial U}{\partial t}+BU = F(U), $$

where

$$B_{i,j} = {\int}_{\Omega}\underline{D}\nabla\phi_{i}\nabla\phi_{j}dx\qquad\forall i,j = 1,2,\dots,Ndof $$

and

$$F(U)_{i} = {\int}_{\Omega}F_{i}(U)\phi_{i}dx\qquad \forall i = 1,2,\dots,Ndof. $$

Ndof is the number of freedom degrees. The sequence (ϕ _i )_{1 ≤ i ≤ Ndof} defines the basis functions of the space V _h.

We use an explicit scheme to discretize time. We denote by N _t the number of time steps. The time step is $dt =\frac {T}{N_{\mathrm {t}}}$. The numerical scheme is:

$$\frac{U^{n}-U^{n-1}}{dt}+B\,U^{n} = F(U^{n-1}),\quad n\in\mathbb{N}^{*} $$

which can be as follows:

$$(I+dt B)U^{n} = U^{n-1}+dt F(U^{n-1}),\quad n\in\mathbb{N}^{*}, $$

where I is the identity matrix and U ⁿ is the solution at time n∗dt.

We give in annex B the Freefem3d code we implement for solving the system.

Appendix B: Algorithm Code Using Freefem3d

The freefem3d code corresponding to the algorithm is the following one:

mesh H = structured((64,64,64),(0,0,0),(1,1,1));

// reading the pore space mesh

function PHI = read(medit,"image.bb",H);

ᅟ

// extract the mesh of the pore space

ᅟ

domain D = domain(PHI>0);

ᅟ

mesh M = tetrahedrize(D,H);

ᅟ

// reaction functions

femfunction f1(M);

femfunction f2(M);

femfunction f3(M);

femfunction f4(M);

femfunction f5(M);

femfunction f6(M);

femfunction ue1(M);

femfunction ue2(M);

femfunction ue3(M);

femfunction ue4(M);

femfunction ue5(M);

femfunction ue6(M);

ᅟ

// parameters double d1 = 1;

double d2 = 1;

double d3 = 1;

double d4 = 1;

double d5 = 1;

double d6 = 1;

double k = 17;

double Kb = 0.0005;

double Km = 0.0005;

double r = 0.2;

double mmu = 1.5;

double nnu = 0;

double zet = 0;

double c1 = 0;

double c2 = 0;

double dt = 0.0416667;

double T = 0;

function Id = 1;

femfunction u1(M) = 6.5⋆10̂-6/int[M](Id);

femfunction u2(M) = 300⋆cos(x⋆y)/int[M](cos(x⋆y)⋆Id);

femfunction u3(M) = 0;

femfunction u4(M) = 0;

femfunction u5(M) = 0;

femfunction u6(M) = 0;

ᅟ

ofstream masse = ofstream("masse.dat");

ᅟ

// save the solutions

masse <<0<<" "<<int[M](u1)<<" "

              <<int[M](u2)<<" " <<int[M](u3)<<" "<<int[M](u4)<<" "

              <<int[M](u5)<<" "<<int[M](u6)<<"\n";

for(double t=1;t<=168;t++)

{

      // save the previous solution

ᅟ

             ue1 = u1;

              ue2 = u2;

              ue3 = u3;

              ue4 = u4;

              ue5 = u5;

              ue6 = u6;

ᅟ

            f1 = k⋆ue1⋆ue2/(Kb+ue1)-(mmu+r+nnu)⋆ue1;

            f2 = -k⋆ue1⋆ue2/(Kb+ue1)+0.5⋆(zet⋆ue5+mmu⋆ue1)

              +(c1⋆ue3+c2⋆ue4)⋆ue5/(Km+ue5);

            f3 = -c1⋆ue3⋆ue5/(Km+ue5)+0.5⋆(zet⋆ue5+mmu⋆ue1);

            f4 = -c2⋆ue4⋆ue5/(Km+ue5);

            f5 = nnu⋆ue1-zet⋆ue5;

            f6 = r⋆ue1;

ᅟ

solve(u1,u2,u3,u4,u5,u6) in M

{

              test(v1,v2,v3,v4,v5,v6)

              int(u1⋆v1) + int(dt⋆d1⋆grad(u1)⋆grad(v1))+

              int(u2⋆v2) + int(dt⋆d2⋆grad(u2)⋆grad(v2))+

              int(u3⋆v3) + int(dt⋆d3⋆grad(u3)⋆grad(v3))+

              int(u4⋆v4) + int(dt⋆d4⋆grad(u4)⋆grad(v4))+

              int(u5⋆v5) + int(dt⋆d5⋆grad(u5)⋆grad(v5))+

              int(u6⋆v6) + int(dt⋆d6⋆grad(u6)⋆grad(v6))

              =

              int((ue1+dt⋆f1)⋆v1)+

              int((ue2+dt⋆f2)⋆v2)+

              int((ue3+dt⋆f3)⋆v3)+

              int((ue4+dt⋆f4)⋆v4)+

              int((ue5+dt⋆f5)⋆v5)+

              int((ue6+dt⋆f6)⋆v6);

       }

      T = T + dt;

ᅟ

       // save the solutions masse <<T<<" "<<int[M](u1)<<" "

              <<int[M](u2)<<" " <<int[M](u3)<<" "<<int[M](u4)<<" "

              <<int[M](u5)<<" "<<int[M](u6)<<"\n"; }

1.1 Appendix C: Reducing the Image Size by Using Octree Method

To compute the solution requires important memory and computing time costs. So, we need to improve calculations cost by reducing the size of the image sample. The method consists in reducing the number of voxels of the image by applying the octree method on the image. In this section, we describe the algorithm written in C++ language. We define the table colourtab such that colourtab[x][y][z] represents the voxels whose coordinates are (x, y, z). dim is the number of subdivisions of each size of the image sample. The statement of colourtab is:

float colourtab[dim][dim][dim];

We define a function CreateCell to implement octree structure from the set of eight voxels. Its inputs are the lower coordinates of the voxels constituting the collection.

void CreateCell(int x,int y,int z)

{

float colour;

ᅟ

// updating data

ᅟ

colour = (tabcolour[x][y][z] + tabcolour[x+1][y][z] +

tabcolour[x+1][y+1][z] + tabcolour[x][y+1][z] +

tabcolour[x][y][z+1] + tabcolour[x+1][y][z+1]+

tabcolour[x+1][y+1][z+1] + tabcolour[x][y+1][z+1])/8;

ᅟ

//saving the results in a file whose identifier is output output

ᅟ

<<x/2<<" "<<y/2<<" "<<z/2<<" "<<(colour>=(1/8)) <<endl;

}

ᅟ

The function CreateOctree creates the approximate image.

void CreateOctree()

 {

    int x, y, z;

    x = 0;

    while(x < dim)

    {

     y = 0;

     while(y < dim)

     {

        z = 0;

        while(z < dim)

        {

                CreateCell(x,y,z);

                z += 2;

        }

        y += 2;

     }

     x += 2;

   }

}