# A Mathematical Model for Transport and Growth of Microbes in Unsaturated Porous Soil

Authors: Ramesh Chandra Timsina; Harihar Khanal; Andrei Ludu; Kedar Nath Uprety
Journal: Mathematical Problems in Engineering
 (2021)

---

## Abstract

In this work, we develop a mathematical model for transport and growth of microbes by natural (rain) water infiltration and flow through unsaturated porous soil along the vertical direction under gravity and capillarity by coupling a system of advection diffusion equations (for concentration of microbes and their growth-limiting substrate) with the Richards equation. The model takes into consideration several major physical, chemical, and biological mechanisms. The resulting coupled system of PDEs together with their boundary conditions is highly nonlinear and complicated to solve analytically. We present both a partial analytic approach towards solving the nonlinear system and finding the main type of dynamics of microbes, and a full-scale numerical simulation. Following the auxiliary equation method for nonlinear reaction-diffusion equations, we obtain a closed form traveling wave solution for the Richards equation. Using the propagating front solution for the pressure head, we reduce the transport equation to an ODE along the moving frame and obtain an analytic solution for the history of bacteria concentration for a specific test case. To solve the system numerically, we employ upwind finite volume method for the transport equations and stabilized explicit Runge–Kutta–Legendre super-time-stepping scheme for the Richards equation. Finally, some numerical simulation results of an infiltration experiment are presented, providing a validation and backup to the analytic partial solutions for the transport and growth of bacteria in the soil, stressing the occurrence of front moving solitons in the nonlinear dynamics.


---

## Body

## 1. Introduction

The flow of water through unsaturated porous soil is of great interest for hydrogeology and water wells monitoring, especially when water flows with contaminated pollutants [1]. With these contaminants, microbial transport in unsaturated porous media is an important aspect from scientific, industrial, and environmental point of view, especially for the water-borne diseases with pathogenic microbes [2]. It is also important for controlling the mass spreading of plant root diseases [3]. The ground water reservoir is recharged from above by infiltrating water coming from various sources [4]. In this stage, water moves downward, from ground surface to the water table, through unsaturated zone. If the surface water contains dissolved contaminants and pollutants, they move downward into the water table by infiltration processes like dispersion and adsorption. These phenomena affect the concentration of pollutants, which move from the surface to the water table [4, 5].Important processes have been identified, and mathematical formulations are developed, enabling scientists to incorporate microbial transport models in flow through porous media [6, 7]. Results of different experiments conducted on soils and other subsurface materials have shown that many environmental factors affect the transport and outcome of microbes in porous media, for example, the chemical composition of interstitial solution, water flow velocity, and physical properties of the solid matrix [6, 8]. Physical processes such as growth, decay, adsorption, straining, advection, motility, diffusion, and dispersion are also known to be involved with the transport and outcome of microbes [9].The transport of microbes in the water phase depends on the flow of water corresponding with advective transport, hydrodynamic dispersion, Brownian diffusion, and motility of some microbes. Microbe movement through a porous medium is also imposed by physical and physicochemical forces, which makes the flow of microbes slow relative to the flow of water. Besides, microbes are subject to decay or death, or their concentration may increase by growth or multiplication, at rates controlled by availability of substances (substrates) that provide nutrients to microbes [9]. Since the consumption of substrate and transport of microbes occur simultaneously, this may also affect the transport phenomenon. To these natural phenomena, we can add human induced mechanisms for retention of microbes, like filtering, adsorption, interception, and sedimentation, the relative importance of which depends on the type of microbes and environmental factors such as ionic strength of solution, toxic chemicals, and properties of the porous medium. By taking into account all these processes, a comprehensive mathematical description for the transport of microbes in unsaturated porous medium can be constructed [10].In this paper, we develop and solve a one-dimensional mathematical model for simulating growth and transport of microbes in unsaturated soil, consisting in a coupled system of four nonlinear partial differential equations. In Section2, we describe the mathematical model for microbial transport and growth through porous medium by coupling a system of convection diffusion equations with Richards equation adopting Darcy’s law. Moreover, the equation for the transport of fluids through unsaturated soils can be derived from the law of mass conservation, by balancing the rate of change of saturation in a closed volume with the rate of change of net flux of fluid into that volume. If the flux of fluid is given by Darcy’s law for flow through a porous medium, the transport equation becomes the Richards equation. There are differences between these two equations. On the one hand, Darcy’s law is valid for any porous medium and any type of pressure drop, while the Richards equation is mainly used for soil fluid dynamics. On the other hand, the hydraulic conductivity parameter K from Richards equation is a function depending on the degree of soil saturation θ (water content in the porous medium) and describes the ease with which the fluid can move through pores under capillary head and gravity. In Darcy’s law, the constant of proportionality is solely a material constant, the porous medium permeability k. If one assumes gravity (g) as the only driving force of the flow, one maps k⟶K=kg/v into the Richards equation, with ν being the kinematic viscosity.We attempt to solve the resulting system of highly nonlinear partial differential equations both analytically and numerically.In Section3, we study the occurrence of front propagating analytic solutions of the main dynamical equation of the model, the Richards nonlinear diffusion equation, which is the only independent equation of the model, in one dependent variable ψ. The other three equations of the model couple all the dynamical variables. Once, the solutions for the volumetric moisture content θz,t and the same thing for the pressure head ψz,t is obtained, one can predict pretty much the dynamics of the rest of the processes and concentrations because the other three dynamical equations of the model functionally depend on ψ. By using several changes of the variable and some linear approximations in the constitutive and material constant equations, we have a weakly dispersive weakly diffusive Burgers type of wave equation for the pressure head ψ and traveling solutions as nontopological solitary waves.For the numerical approximation, we employ finite difference method for the Richards equation and finite volume method for the transport equations. In our earlier work [11], we studied performance of various explicit and implicit finite difference schemes for the Richards equation. A numerical procedure based on upwind finite volume discretization of the transport equations and their boundary conditions is developed in Section 4. The closed form analytic solution obtained in Section 3 contains several parameters (constants of integration). Variations in these parameters correspond to variations in the physical constants (material properties) present in the mathematical model, thereby representing change in different experimental setup. By hits and trial approach, we obtained one traveling wave solution with a specific wave front velocity which mimics the qualitative behavior of the solution of the nonlinear system. The authors believe that this soliton-like solution for the transport model raises an interesting problem for the researchers to look at the solution using data science/machine learning approach and try to find the optimal values for the constants in the analytic solution that represents a given set of physical parameters of the model. Some numerical simulation results are presented in Section 5. Finally in Section 6, we conclude our results.


## 2. The Mathematical Model

### 2.1. Flow in Unsaturated Porous Medium

Flow in unsaturated soils is a special type of flow in porous media when the void spaces are not completely filled with fluid. Hydraulic processes at land surface and subsurface, such as infiltration, drainage, evaporation, capillary elevation of deep-level water, and absorption of soil-water by plant roots zone, all can be reduced to unsaturated flow problems [4]. Virtually, in all studies performed for unsaturated regime, the flow is assumed to obey the classical Richards equation [12], which has been directly obtained from applying the equation of continuity to the Darcy-Buckingham flow law [13]. The Richards equation in one-dimensional vertical z− direction (z>0 downward) with no source/sink term has the form(1)∂θ∂t=∂∂zK∂ψ∂z−∂K∂z,where θz,t is the volumetric moisture content, ψz,t is the pressure head, and Kψ is the unsaturated hydraulic conductivity. To solve (1), we use two empirical constitutive relations, one for the moisture content and the other for the hydraulic conductivity [14]:(2)θψ=θr+αθs−θrα+ψβ,Kψ=KsAA+ψγ,where θs and θr represent the saturated and residual moisture content, respectively; Ks corresponds to the saturated hydraulic conductivity; and A, α, β, and γ are dimensionless soil parameters.


### 2.2. Microbial Transport in the Porous Medium

The mathematical equation for transport of microbes in unsaturated porous media is based on the continuity or mass balance equation. In an elementary volume of porous medium, it states that the rate of change of the total biomass equals the divergence of flux plus the rate of contribution to the biomass (growth) minus the rate of removal of the biomass (decay) [9].LetMz,t denote the concentration of microbes in the water phase defined as the mass of microbial cells per unit volume of water, ρs the bulk density of the porous medium, and Maz,t the mass of microbes attached to the solid phase. At moment t, the total biomass in the elementary volume is the sum of the mass in the water phase θM plus the mass in the attached phase ρsMa [6, 8]. The governing equation for microbial transport in the porous media in one dimension is given by the following partial differential equation:(3)∂θM+ρsMa∂t=∂∂z−vwθM−vθM+τθD∂M∂z+Rg−Rd,where vw is the velocity of water flow, v is the chemostatic velocity of the microbes, τ is the tortuosity of the pore space, D is a diffusion coefficient, Rg is the growth rate of microbes combined for water and attached biomass phase (mass of microbes per unit volume and time), and Rd is the decay rate of the microbes combined for water and attached biomass phase.


### 2.3. Transport Equation for Growth-Limiting Substrate

In a porous medium, the transport of growth-limiting substrate and transport of microbes occur simultaneously. Therefore, the growth of microbes depends upon the flow of substrate. Moreover, the microbe transport and substrate transport also occur at the same time, so there is a continuous change of substrate concentration [9]. Thus, it is necessary to derive and solve the substrate transport equation with that for microbial transport to obtain the substrate concentration. Processes such as advection with the flowing water, diffusion and dispersion phenomena relative to water in response to concentration gradients, sorption by specific adsorption reactions, and consumption of substrate by microbes affect the transport of substrate in porous media [6, 8]. Including all these, the transport equation in one dimension for biologically reactive substrate can be written as(4)∂θS+ρsSa∂t=∂∂z−θvwS+θDs∂S∂z−Rb,where Sz,t is the mass of substrate per unit volume of water, Saz,t is the mass of substrate adsorbed per unit mass of solid matrix, Ds is the diffusion-dispersion coefficient of the substrate, and Rb is the rate at which the microbes consume substrate combined for the water and attached biomass phases [9].


### 2.4. Microbial Behavior in Porous Media

We use the Monod model [15] for microbial growth and substrate consumption. Let Rgw, Rga, Rdw, and Rda denote the growth and decay rates of microbes in water and attached biomass, respectively. Then, the growth rate of microbes combined for water and attached biomass phase (mass of microbes per unit volume and time) is defined as(5)Rg=Rgw+Rga=μθM+μρsMa,where μ is the specific growth rate of microbes given in [15] in the form(6)μ=μmSks+S,where μm is specific growth rate of microbes and ks is saturation constant. The decay rate of microbes is defined as(7)Rd=Rdw+Rda=bθM+bρsMa,where b is the decay or inactivation constant. The substrate consumption Rb is given by(8)Rb=μYθM+ρsMa,where Y is defined as the maximum yield coefficient, that is, the amount of biomass produced per unit mass of substrate consumed.Various absorption models can be used to eliminate one of the dependent variablesM,Ma in (2.4) and S,Sa in (4). For the substrate in the attached biomass, we assume equilibrium adsorption isotherm (Freundlich) [16] as(9)Sa=ρskfsSNs,where kfs is Freundlich isotherm slope parameter. For the microbial growth in the attached biomass phase, we use the following nonequilibrium microbial kinetics:(10)∂ρsMa∂t=Ra−Ry+Rga−Rda,where Ra is the rate of microbial attachment and Ry is the rate of microbial detachment. Following [6], we define(11)Ra=ka1−MaMmaθM,Ry=kyρsMa,where Mma is the maximum concentration of microbes in attached phase, ka is attachment coefficient, and ky is detachment coefficient.


### 2.5. Mathematical Model

By takingNs=1 in (9) and coupling Richards (1) for porous flow, with (4) for the transport of biologically reactive solutes, with (10) for the microbial kinetics, and with (3) for microbial transport in porous media, we can write our model as a 4×4 nonlinear PDE system in the functions Ψ,M,Ma, and S in a bounded region 0<z<L as follows:(12)∂θ∂t=∂∂zKψ∂ψ∂z−∂K∂z,(13)∂θ+ρskfsS∂t=∂∂z−vwθS+θDs∂S∂z−μmSks+SYθM+ρsMa,(14)∂Ma∂t=kaMma−MaρsMmaθM−kyMa+μmSks+S−bMa,(15)∂θM∂t=∂∂z−vwθM−vθM+τθD∂M∂z+kyρsMa−ka1−MaMmaθM+μmSks+S−bθM,subject to the following initial and boundary conditions:(i)
Initial condition:ψz,0=ψ0z, Mz,0=M0z, Maz,0=Maz, Sz,0=S0z.(ii)
Lower boundary condition:ψL,t=βt, ∂M/∂zL,t=0, ∂S/∂zL,t=0.(iii)
Upper boundary condition:(16)∂ψ∂z0,t=qt,if0<t<tα,0,ift>tα,where the functions of the initial condition and boundary conditions are chosen correspondingly.


## 2.1. Flow in Unsaturated Porous Medium

Flow in unsaturated soils is a special type of flow in porous media when the void spaces are not completely filled with fluid. Hydraulic processes at land surface and subsurface, such as infiltration, drainage, evaporation, capillary elevation of deep-level water, and absorption of soil-water by plant roots zone, all can be reduced to unsaturated flow problems [4]. Virtually, in all studies performed for unsaturated regime, the flow is assumed to obey the classical Richards equation [12], which has been directly obtained from applying the equation of continuity to the Darcy-Buckingham flow law [13]. The Richards equation in one-dimensional vertical z− direction (z>0 downward) with no source/sink term has the form(1)∂θ∂t=∂∂zK∂ψ∂z−∂K∂z,where θz,t is the volumetric moisture content, ψz,t is the pressure head, and Kψ is the unsaturated hydraulic conductivity. To solve (1), we use two empirical constitutive relations, one for the moisture content and the other for the hydraulic conductivity [14]:(2)θψ=θr+αθs−θrα+ψβ,Kψ=KsAA+ψγ,where θs and θr represent the saturated and residual moisture content, respectively; Ks corresponds to the saturated hydraulic conductivity; and A, α, β, and γ are dimensionless soil parameters.


## 2.2. Microbial Transport in the Porous Medium

The mathematical equation for transport of microbes in unsaturated porous media is based on the continuity or mass balance equation. In an elementary volume of porous medium, it states that the rate of change of the total biomass equals the divergence of flux plus the rate of contribution to the biomass (growth) minus the rate of removal of the biomass (decay) [9].LetMz,t denote the concentration of microbes in the water phase defined as the mass of microbial cells per unit volume of water, ρs the bulk density of the porous medium, and Maz,t the mass of microbes attached to the solid phase. At moment t, the total biomass in the elementary volume is the sum of the mass in the water phase θM plus the mass in the attached phase ρsMa [6, 8]. The governing equation for microbial transport in the porous media in one dimension is given by the following partial differential equation:(3)∂θM+ρsMa∂t=∂∂z−vwθM−vθM+τθD∂M∂z+Rg−Rd,where vw is the velocity of water flow, v is the chemostatic velocity of the microbes, τ is the tortuosity of the pore space, D is a diffusion coefficient, Rg is the growth rate of microbes combined for water and attached biomass phase (mass of microbes per unit volume and time), and Rd is the decay rate of the microbes combined for water and attached biomass phase.


## 2.3. Transport Equation for Growth-Limiting Substrate

In a porous medium, the transport of growth-limiting substrate and transport of microbes occur simultaneously. Therefore, the growth of microbes depends upon the flow of substrate. Moreover, the microbe transport and substrate transport also occur at the same time, so there is a continuous change of substrate concentration [9]. Thus, it is necessary to derive and solve the substrate transport equation with that for microbial transport to obtain the substrate concentration. Processes such as advection with the flowing water, diffusion and dispersion phenomena relative to water in response to concentration gradients, sorption by specific adsorption reactions, and consumption of substrate by microbes affect the transport of substrate in porous media [6, 8]. Including all these, the transport equation in one dimension for biologically reactive substrate can be written as(4)∂θS+ρsSa∂t=∂∂z−θvwS+θDs∂S∂z−Rb,where Sz,t is the mass of substrate per unit volume of water, Saz,t is the mass of substrate adsorbed per unit mass of solid matrix, Ds is the diffusion-dispersion coefficient of the substrate, and Rb is the rate at which the microbes consume substrate combined for the water and attached biomass phases [9].


## 2.4. Microbial Behavior in Porous Media

We use the Monod model [15] for microbial growth and substrate consumption. Let Rgw, Rga, Rdw, and Rda denote the growth and decay rates of microbes in water and attached biomass, respectively. Then, the growth rate of microbes combined for water and attached biomass phase (mass of microbes per unit volume and time) is defined as(5)Rg=Rgw+Rga=μθM+μρsMa,where μ is the specific growth rate of microbes given in [15] in the form(6)μ=μmSks+S,where μm is specific growth rate of microbes and ks is saturation constant. The decay rate of microbes is defined as(7)Rd=Rdw+Rda=bθM+bρsMa,where b is the decay or inactivation constant. The substrate consumption Rb is given by(8)Rb=μYθM+ρsMa,where Y is defined as the maximum yield coefficient, that is, the amount of biomass produced per unit mass of substrate consumed.Various absorption models can be used to eliminate one of the dependent variablesM,Ma in (2.4) and S,Sa in (4). For the substrate in the attached biomass, we assume equilibrium adsorption isotherm (Freundlich) [16] as(9)Sa=ρskfsSNs,where kfs is Freundlich isotherm slope parameter. For the microbial growth in the attached biomass phase, we use the following nonequilibrium microbial kinetics:(10)∂ρsMa∂t=Ra−Ry+Rga−Rda,where Ra is the rate of microbial attachment and Ry is the rate of microbial detachment. Following [6], we define(11)Ra=ka1−MaMmaθM,Ry=kyρsMa,where Mma is the maximum concentration of microbes in attached phase, ka is attachment coefficient, and ky is detachment coefficient.


## 2.5. Mathematical Model

By takingNs=1 in (9) and coupling Richards (1) for porous flow, with (4) for the transport of biologically reactive solutes, with (10) for the microbial kinetics, and with (3) for microbial transport in porous media, we can write our model as a 4×4 nonlinear PDE system in the functions Ψ,M,Ma, and S in a bounded region 0<z<L as follows:(12)∂θ∂t=∂∂zKψ∂ψ∂z−∂K∂z,(13)∂θ+ρskfsS∂t=∂∂z−vwθS+θDs∂S∂z−μmSks+SYθM+ρsMa,(14)∂Ma∂t=kaMma−MaρsMmaθM−kyMa+μmSks+S−bMa,(15)∂θM∂t=∂∂z−vwθM−vθM+τθD∂M∂z+kyρsMa−ka1−MaMmaθM+μmSks+S−bθM,subject to the following initial and boundary conditions:(i)
Initial condition:ψz,0=ψ0z, Mz,0=M0z, Maz,0=Maz, Sz,0=S0z.(ii)
Lower boundary condition:ψL,t=βt, ∂M/∂zL,t=0, ∂S/∂zL,t=0.(iii)
Upper boundary condition:(16)∂ψ∂z0,t=qt,if0<t<tα,0,ift>tα,where the functions of the initial condition and boundary conditions are chosen correspondingly.


## 3. Traveling Front Solutions for the Richards Equation

The nonlinear system equations (12)–(15) have a particular feature which allows us to introduce some analytical considerations. Richards’ equation (12) is independently decoupled for the other three equations, and once we use a certain porous fluid material model, for example, the constitutive equations (2), (12) becomes a nonlinear PDE in one function ψz,t. Once a solution is found for this equation, by implementing it into the remaining three equations, (13)–(15), they become a system of linear PDE with variable coefficients for which a variety of solutions can be obtained.Consequently, in this section, we study some properties of the Richards equation. By dividing (12) by ∂θ/∂ψ, we can rewrite this equation in the form(17)∂ψ∂t=θ1K1∂ψ∂z2+θ1K∂2ψ∂z2−θ1K1∂ψ∂z,where we consider K=Kψ and θ=θψ to be invertible differentiable functions, at least for conveniently chosen intervals for the K,ψ,θ variables, and where we denote(18)θ1ψ=∂ψ∂θ−1,K1ψ=∂K∂ψ.Equation (17) is a Hamilton–Jacobi type of equation, and under corresponding Cauchy initial conditions and appropriate boundary conditions, it has a unique solution [17]. To understand more about such a solution, we notice that for slow variation coefficient functions, like the linear approximations θ≃C1+C2ψ and K≃C3+C4ψ, (17) becomes a Burgers equation:(19)ηt+C2C4ηz=3C2C4ηηz+C2C3ηzz+C2C4ηηzz,in η=ψz, where suffixes represent derivations. This result arrives somehow in a natural way because (17) is a Richards (or Darcy) equation which is a particular case of the Navier–Stokes equation, and in the limit of a porous homogeneous medium, the Navier–Stokes equation turns into some type of nonlinear evolution equation like Burgers, Korteweg-de Vries, or generalizations of these depending on the boundary conditions. Usually, the analytical solutions of the Burgers equations evolve into shock formation and blow-up solutions. In the case of the Richards equation, the sign in front of the time derivative is opposite to the sign of the corresponding term in the Burgers equation, so the front solutions tend to stabilize in time.In the following, we introduce a transformationψ=ψσ fulfilling the condition(20)dψdσ=σ0Kσ,where we assume the function Kψ=Kψσ is invertible, differentiable, and nonzero and where σ0 is an arbitrary constant. With this type of mapping ψ⟶σ, (17) becomes(21)∂σ∂t=θ1K∂2σ∂z2−θ1K1dψdσ−1∂σ∂z,which is a diffusion equation with nonlinearity resulting from its functional coefficients. For example, if K is given by the material relation in (3), by integration of (20), we have(22)σ=σ1+KsAψσ0,F211,1γ,1+1γ,−ψγ,which is a strictly monotonic function. For small γ>0, (22) generates a linear dependence σ∼ψ, and for larger γ⟶∞, (22) approaches a tanh function.An exact solution for the Richards equation can be obtained in the case of slow variation of the functional coefficients with respect to the new variableσ. By using such a linear approximation(23)θ1K≃c1+c2σ,θ1K1dψdσ−1≃c3+c4σ,we obtain for (21) an exact traveling front solution [18]:(24)z−Vt=z0+c22c3−V∫1σu0+c2c4c3−Vu+c2lnc1+c2uV−c3+c1c4V−c3lnc1+c2uV−c3−1du,where ξ=z−Vt is the comoving coordinate; V is a free parameter representing the velocity of the front; and z0,u0 are constants of integration. For example, if we choose a simple case where β=γ=1, the following results:(25)σ=σ1+KsAψσ0.Consequently, in the first order of their Taylor series, we have(26)θ1K1=−AKsαθr−θsA+σ02α+σ02+2A2Ks2αθr−θsA+α+2σ0σ0A+σ03α+σ03σ+O2σ,θ1K1dψdσ−1=−αθr−θsσ0A+σ02α+σ02+2AKsαθr−θsA+α+2σ0A+σ03α+σ03σ+O2σ.In continuation, we can choose the arbitrary parametersσ0,σ1 so that (24) will match appropriate boundary and initial conditions. In particular, we can always choose the change of variable ψσ such that the moisture θ has an extremum at σ=0 which would represent the choice of a certain value for the pressure head ψ0 generating maximum, or minimum, moisture content. Consequently, θ10K10=0 which means a choice c3=0. In this case, the solution equation (24) of (21) has an even simple form:(27)z−Vt=z0−ec2σ0Ei−c2σ0+lnc1+c2σ,where σ0 is constant of integration and Ei is the exponential integral function. Because the inverse of the exponential integral function is multivalued, the solutions σz−Vt from (27) can be classified into two types of traveling front solutions: wetting fronts and dewatering fronts.Equation (17) is a typical diffusion-reaction equation [19], and it is integrable with soliton solutions. If we consider a traveling front solution propagating with velocity V, (17) can be written in the comoving frame ζ=z−Vt in the form(28)dψdζ+Ad2ψdζ2+Bdψdζ2=0,with(29)Aψ=θ1KV−θ1K1,Bψ=θ1K1V−θ1K1.We can write the differential equation fulfilled by the inverse functionζ=ζψ in the form(30)d2ζdψ2−1Adζdψ+BA=0,which has solution(31)ζψ=C1+∫0ψζ1sC2−∫0sBpApζ1pdpds,where(32)ζ1s=exp∫0sdpAp,and the traveling solution ψz−Vt is obtained by inversion of (30).In the following, we plug this propagating front solution for the pressure head in (15) to obtain the evolution of the concentration of microbes in water phase Mz,t. We consider in a first approximation that Ma≃0, and we take a realistic value for γ=4.75 (see Section 5) and for the other constants of material. In the comoving frame ζ=z−Vt, (15) becomes a linear ODE with variable coefficients.


## 4. Numerical Methods

### 4.1. Numerical Method for the Richards Equation

For the numerical approximation to the solution of (12), following [11, 20], we use the following equivalent Kirchhoff transformed equation:(33)cϕ∂ϕ∂t=∂2ϕ∂z2,where the functional coefficient c is given by(34)cϕ=αβθs−θrhβ−1K¯hα+|h|β2,which depends on the pressure head ψ through the Kirchhoff transformation defined as(35)ϕh=∫0hK¯λdλ,with h=ψ−z and K¯h=Kψ. We solve the above one-dimensional Kirchhoff transformed Richards equation (33) numerically using a stabilized explicit Runge–Kutta–Legendre super-time-stepping (RKL) finite difference scheme. Detailed description of the numerical procedure and stability analysis can be found in our earlier work [11]. The numerical method is unconditionally stable and easily parallelizable and can be implemented in 2D and 3D geometry.


### 4.2. Numerical Method for Microbial Transport and Growth Equations

Now, we discretize the transport equations (13) and (15) and the growth equation (14) using the finite volume method [21, 22].


#### 4.2.1. Finite Volume Discretization

We partition the computational domain (the slab of lengthL) into Nz finite cells Zj=zj−1/2,zj+1/2,j=1,…,Nz, with length Δz=L/Nz and define the nodes zj as the midpoint of the cell Zj. That is, the nodes zj and the faces zj±1/2 of the control volumes Zj are given by(36)z0=z1−12=0,zj−12=jΔz,zj=12zj−12+zj+12,j=1,…,Nz,zNz+1=zNz+12=L.LetΔt=T/Nt be the time increments and define the discrete time-steps tn=nΔt, n=0,1,2,…,Nt.The initial time is denoted byn=0, and the discrete approximation for S at the grid zj,tn is denoted by Sjn. We also regard Sjn as an approximation to mean value of Sz,t in the cell Zj at time tn as(37)Sjn=1Δz∫zj−1/2zj+1/2Sz,tndz.The variablesMjn and θjn can be derived in similar manner (13)If an approximate solutionSjn is assumed to be known at all grid points at time tn, a method must be specified to advance the solution to time tn+1, subject to the boundary conditions. For this, in finite volume method, we integrate equation (14) into the finite cell Zj and approximate the integral to get the numerical scheme.(38)∫zj−1/2zj+1/2∂∂tθ+ρskfsS−∂∂z−vwθS+θDs∂S∂z+μmSks+SYθM+ρsMadz=0.The first term in the integral (38) is approximated using forward difference in time as(39)∫zj−1/2zj+1/2∂θ+ρskfS∂tdz≈Δzθjn+1+ρskfSjn+1−θjn+ρskfSjnΔt.By applying the Gauss divergence theorem, the second term in the integral (38) can be reduced to sum fluxes at the boundary faces of the finite volume Zj. Clearly, these fluxes do have two parts, advection and diffusion. For the advective fluxes we use upwind approximation.(40)∫zj−1/2zj+1/2∂∂z−vwθSdz≈−vwθSj+1/2n−θSj−1/2n≈−vwθSj+1n−θSjn2−θSj−1n−θSjn2,∫zj−1/2zj+1/2∂∂zθDs∂S∂zdz≈Dsθ∂S∂zj+1/2n−θ∂S∂zj−1/2n≈Dsθj+1n+θjn2Sj+1n+SjnΔz−θj−1n+θjn2Sjn+Sj−1nΔz.Approximating the third term as(41)∫zj−1/2zj+1/2μmSYks+SθM+ρsMadz≈ΔzμmSjnθjnMjn+ρsMjanYks+Sjn,yields the following scheme:(42)Sjn+1=1θjn+1+ρskfθjn+ρskfSjn+ΔtΔz−vwθSj+1n−θSjn2−θSj−1n−θSjn2+Dsθj+1n+θjn2Sj+1n+SjnΔz−θj−1n+θjn2Sjn+Sj−1nΔz+ΔzμmSjnθjnMjn+ρsMjanYks+Sjn.Similarly, we write the finite volume discretization for (14) and (15).Discretization of the growth model equation (14) is as follows:(43)Mjan+1=Mjan+Δtkaρs1−MjanMmaθMjn−kyMjan+μmSjnks+Sjn−bMjan.Discretization of the microbial transport equation (15) is as follows:(44)Mjn+1=1θjn+1θjnMjn+ΔtΔz−vw−vxθMj+1n−θMjn2−θMj−1n−θMjn2+τDθj+1n+θjn2Mj+1n−MjnΔz−θj−1n+θjn2Mjn−Mj−1nΔz+ΔtkyρsMjan+1−ka1−Mjan+1MmaθjnMjn+μmSjnks+SjnθjnMjn.


## 4.1. Numerical Method for the Richards Equation

For the numerical approximation to the solution of (12), following [11, 20], we use the following equivalent Kirchhoff transformed equation:(33)cϕ∂ϕ∂t=∂2ϕ∂z2,where the functional coefficient c is given by(34)cϕ=αβθs−θrhβ−1K¯hα+|h|β2,which depends on the pressure head ψ through the Kirchhoff transformation defined as(35)ϕh=∫0hK¯λdλ,with h=ψ−z and K¯h=Kψ. We solve the above one-dimensional Kirchhoff transformed Richards equation (33) numerically using a stabilized explicit Runge–Kutta–Legendre super-time-stepping (RKL) finite difference scheme. Detailed description of the numerical procedure and stability analysis can be found in our earlier work [11]. The numerical method is unconditionally stable and easily parallelizable and can be implemented in 2D and 3D geometry.


## 4.2. Numerical Method for Microbial Transport and Growth Equations

Now, we discretize the transport equations (13) and (15) and the growth equation (14) using the finite volume method [21, 22].


### 4.2.1. Finite Volume Discretization

We partition the computational domain (the slab of lengthL) into Nz finite cells Zj=zj−1/2,zj+1/2,j=1,…,Nz, with length Δz=L/Nz and define the nodes zj as the midpoint of the cell Zj. That is, the nodes zj and the faces zj±1/2 of the control volumes Zj are given by(36)z0=z1−12=0,zj−12=jΔz,zj=12zj−12+zj+12,j=1,…,Nz,zNz+1=zNz+12=L.LetΔt=T/Nt be the time increments and define the discrete time-steps tn=nΔt, n=0,1,2,…,Nt.The initial time is denoted byn=0, and the discrete approximation for S at the grid zj,tn is denoted by Sjn. We also regard Sjn as an approximation to mean value of Sz,t in the cell Zj at time tn as(37)Sjn=1Δz∫zj−1/2zj+1/2Sz,tndz.The variablesMjn and θjn can be derived in similar manner (13)If an approximate solutionSjn is assumed to be known at all grid points at time tn, a method must be specified to advance the solution to time tn+1, subject to the boundary conditions. For this, in finite volume method, we integrate equation (14) into the finite cell Zj and approximate the integral to get the numerical scheme.(38)∫zj−1/2zj+1/2∂∂tθ+ρskfsS−∂∂z−vwθS+θDs∂S∂z+μmSks+SYθM+ρsMadz=0.The first term in the integral (38) is approximated using forward difference in time as(39)∫zj−1/2zj+1/2∂θ+ρskfS∂tdz≈Δzθjn+1+ρskfSjn+1−θjn+ρskfSjnΔt.By applying the Gauss divergence theorem, the second term in the integral (38) can be reduced to sum fluxes at the boundary faces of the finite volume Zj. Clearly, these fluxes do have two parts, advection and diffusion. For the advective fluxes we use upwind approximation.(40)∫zj−1/2zj+1/2∂∂z−vwθSdz≈−vwθSj+1/2n−θSj−1/2n≈−vwθSj+1n−θSjn2−θSj−1n−θSjn2,∫zj−1/2zj+1/2∂∂zθDs∂S∂zdz≈Dsθ∂S∂zj+1/2n−θ∂S∂zj−1/2n≈Dsθj+1n+θjn2Sj+1n+SjnΔz−θj−1n+θjn2Sjn+Sj−1nΔz.Approximating the third term as(41)∫zj−1/2zj+1/2μmSYks+SθM+ρsMadz≈ΔzμmSjnθjnMjn+ρsMjanYks+Sjn,yields the following scheme:(42)Sjn+1=1θjn+1+ρskfθjn+ρskfSjn+ΔtΔz−vwθSj+1n−θSjn2−θSj−1n−θSjn2+Dsθj+1n+θjn2Sj+1n+SjnΔz−θj−1n+θjn2Sjn+Sj−1nΔz+ΔzμmSjnθjnMjn+ρsMjanYks+Sjn.Similarly, we write the finite volume discretization for (14) and (15).Discretization of the growth model equation (14) is as follows:(43)Mjan+1=Mjan+Δtkaρs1−MjanMmaθMjn−kyMjan+μmSjnks+Sjn−bMjan.Discretization of the microbial transport equation (15) is as follows:(44)Mjn+1=1θjn+1θjnMjn+ΔtΔz−vw−vxθMj+1n−θMjn2−θMj−1n−θMjn2+τDθj+1n+θjn2Mj+1n−MjnΔz−θj−1n+θjn2Mjn−Mj−1nΔz+ΔtkyρsMjan+1−ka1−Mjan+1MmaθjnMjn+μmSjnks+SjnθjnMjn.


## 4.2.1. Finite Volume Discretization

We partition the computational domain (the slab of lengthL) into Nz finite cells Zj=zj−1/2,zj+1/2,j=1,…,Nz, with length Δz=L/Nz and define the nodes zj as the midpoint of the cell Zj. That is, the nodes zj and the faces zj±1/2 of the control volumes Zj are given by(36)z0=z1−12=0,zj−12=jΔz,zj=12zj−12+zj+12,j=1,…,Nz,zNz+1=zNz+12=L.LetΔt=T/Nt be the time increments and define the discrete time-steps tn=nΔt, n=0,1,2,…,Nt.The initial time is denoted byn=0, and the discrete approximation for S at the grid zj,tn is denoted by Sjn. We also regard Sjn as an approximation to mean value of Sz,t in the cell Zj at time tn as(37)Sjn=1Δz∫zj−1/2zj+1/2Sz,tndz.The variablesMjn and θjn can be derived in similar manner (13)If an approximate solutionSjn is assumed to be known at all grid points at time tn, a method must be specified to advance the solution to time tn+1, subject to the boundary conditions. For this, in finite volume method, we integrate equation (14) into the finite cell Zj and approximate the integral to get the numerical scheme.(38)∫zj−1/2zj+1/2∂∂tθ+ρskfsS−∂∂z−vwθS+θDs∂S∂z+μmSks+SYθM+ρsMadz=0.The first term in the integral (38) is approximated using forward difference in time as(39)∫zj−1/2zj+1/2∂θ+ρskfS∂tdz≈Δzθjn+1+ρskfSjn+1−θjn+ρskfSjnΔt.By applying the Gauss divergence theorem, the second term in the integral (38) can be reduced to sum fluxes at the boundary faces of the finite volume Zj. Clearly, these fluxes do have two parts, advection and diffusion. For the advective fluxes we use upwind approximation.(40)∫zj−1/2zj+1/2∂∂z−vwθSdz≈−vwθSj+1/2n−θSj−1/2n≈−vwθSj+1n−θSjn2−θSj−1n−θSjn2,∫zj−1/2zj+1/2∂∂zθDs∂S∂zdz≈Dsθ∂S∂zj+1/2n−θ∂S∂zj−1/2n≈Dsθj+1n+θjn2Sj+1n+SjnΔz−θj−1n+θjn2Sjn+Sj−1nΔz.Approximating the third term as(41)∫zj−1/2zj+1/2μmSYks+SθM+ρsMadz≈ΔzμmSjnθjnMjn+ρsMjanYks+Sjn,yields the following scheme:(42)Sjn+1=1θjn+1+ρskfθjn+ρskfSjn+ΔtΔz−vwθSj+1n−θSjn2−θSj−1n−θSjn2+Dsθj+1n+θjn2Sj+1n+SjnΔz−θj−1n+θjn2Sjn+Sj−1nΔz+ΔzμmSjnθjnMjn+ρsMjanYks+Sjn.Similarly, we write the finite volume discretization for (14) and (15).Discretization of the growth model equation (14) is as follows:(43)Mjan+1=Mjan+Δtkaρs1−MjanMmaθMjn−kyMjan+μmSjnks+Sjn−bMjan.Discretization of the microbial transport equation (15) is as follows:(44)Mjn+1=1θjn+1θjnMjn+ΔtΔz−vw−vxθMj+1n−θMjn2−θMj−1n−θMjn2+τDθj+1n+θjn2Mj+1n−MjnΔz−θj−1n+θjn2Mjn−Mj−1nΔz+ΔtkyρsMjan+1−ka1−Mjan+1MmaθjnMjn+μmSjnks+SjnθjnMjn.


## 5. Numerical Results

### 5.1. Simulation Setup

Richards equation for unsaturated flow was linked to the microbial transport equation along with transport of growth-limiting substrate. In this simulation, we consider a vertical soil column of depthL=70 cm in a time period of tmax=1hr hr. The simulation starts with a uniform saturation θ=0.1 cm 3/cm3, and a constant water head ψ=−61.5 cm is maintained at the bottom boundary z=L. At the upper boundary z=0 (the soil surface), a constant flux qt=13.69 cm/hr for t<0.3 hr and zero normal flux condition for t>0.3 hr is adopted. Equal amounts of 103 cfu/ml microbial concentration in water phase and attached to solid matrix with 10 mg/100 ml of substrate concentration were set through column at 0.34 cm/sec of water flow velocity. Since no experimental evidence of the model is reported, all fixed model parameter values for soil characteristic used for model simulation are taken from Haverkamp et al. (1977). The values for the empirical parameters θs, θr, ks, ψ for soil type are adopted from [14]. In addition, the parameters related to the microbial transport and fate are listed in Table 1 ([6, 8]).Table 1
Kinetic parameters.
ParameterValueDescriptionD1.51×10−10m2/sDiffusion coefficientτ2.0/3.0Tortuosityvw0.34 cm/sSpeed of water flowv0.66 cm/sSpeed of microbe flowμm0.94 microbes/sSpecific growth rate of microbesks100 cm/sSaturation constantYs6.5 g/mol ATPBiomass produced per unit mass of substrate consumedB1.5×10−7microbes/hrDeath rate of bacteriaMma1100Maximum concentration of microorganisms in attached phaseka0.85Attachment coefficientkf0.95Freundlich isotherm slope parameterky0.37Detachment coefficientρs1Density of solid matrixThe moisture content has been taken as a function of the pressure head (suction); the constitutive relationship betweenθz,t and ψz,t is(45)θψ=θr+αθs−θrα+ψβ=0.075+1.611×106.287−0.0751.611×106+ψ3.96,Kψ=KsAA+ψγ=34×1.175×1061.175×106+ψ4.74.The simulation results are illustrated in Figures1–7. We show the variation trend of microbial concentration suspended in water Mz,t and attached to solid matrix Maz,t and the growth-limiting substrate Sz,t in depth.Figure 1
Front propagation analytic solution for the time evolution of the concentration of microbes in water phaseMz,t, plotted in relative units for several equidistantly separated depths z. The solution is obtained by solving (15) for M with θ given by the solution in (31), Ma approximated zero, γ=4.75, and the constants of material conveniently chosen. In the inset, we present the corresponding solutions σψ of (27) for several values for γ for comparison.Figure 2
Microbe concentration in waterMz,t (15) plotted vs. depth z for four moments of time. It can be seen from figure that the concentration of microbes renders due to the flow process and the consumption of growth limiting substrate but with decreasing in moisture content and because the property of aw microbial growth in stationary phase.Figure 3
Substrate concentration in waterSz,t plotted vs. depth z for four moments of time. It can be seen that the concentration of growth-limiting substrate decreases rapidly because of high consumption by microbes and increases apparently due to the decrease in moisture of the soil.Figure 4
Graphical representation of transport of microbes in aqueous phase for different depths (a) and substrates (b). Different colored curves represent distributions taken at moments of time separated by 0.1 hrs.
(a)(b)Figure 5
Graphical representation of transport of microbes in aqueous phase for different time periods (a), substrates, (b) and microbes in attached phase (c).
(a)(b)(c)Figure 6
Variational trend of microbes in aqueous phase in depth and time.Figure 7
Variational trend of substrate in aqueous phase in depth and time.The numerical procedure developed in the previous section is written in Python 3.7 Compiler with Intel® Core™ i5-2450M CPU @ 2.50 GHz.


### 5.2. Results and Discussion

The proposed model was applied to simulate the transport and growth of microbes with the transport of growth-limiting substrate through an unsaturated finite column. The growth of microorganisms is often controlled by a single substrate and the specific growth rate is assumed to be a function of the concentration of that rate limiting substrate. It also depends on the behavior of moisture, i.e., wetting and drainage, and the available substrate property, i.e., water activityaw of substrate. The aw of substrate is the ratio of vapor pressure of the substrate to the vapor pressure of pure water. Decrease in moisture content leads to decrease in aw value of the substrate that either decreases or ceases the microbial growth. The movement of microorganisms depends on the flow velocity, dispersion, diffusion, Brownian motion, and also chemostatic and tumbling motion of microbes. Here, in our mathematical model, we used the hypothetical values for microbes and growth-limiting substrate at the top of the soil column, and then the total microbial count, moisture content count, and substrate content count down the soil column were calculated. The model shows that the initial bacterial count, i.e., 1.1×103 cfu/ml, increases down the soil column in presence of substrate, i.e., 10 mg/100 ml, with the moisture of the porous medium. However, as the microbes and nutrients percolate down the column, there is steady decrease in moisture. The decrease in moisture showed direct influence on microbial growth, and the total microbial count decreased overall. Likewise, the substrate concentration showed decrease in concentration at the upper part of the column; however, the substrate appears to apparently increase towards the end due to reduction of moisture in soil. The decrease in microbial count down the column by continuous decrease in moisture results in adverse condition for microbial growth. Likewise, the apparent increase in substrate down the column is caused by decrease in moisture that concentrates the substrate in a limited amount of moisture. These phenomena are interpreted in the figures.Figure8 shows the front propagation analytic solution for the time evolution of the concentration of microbes in water phase. Figure 1 explains the growth of microbes in aqueous phase with respect to depth for flow velocity, whereas Figure 2 explains the growth of microorganisms attached to solid matrix with respect to depth for flow velocity. Furthermore, Figure 3 shows the consumption behavior of substrate. Figure 4 explains the transport of microbes and substrate for different time in aqueous and attached phase, and Figures 5–7 describe the variational trend of microbes and substrate in aqueous phase and attached phase in depth and time. The continuous decrease in moisture thus decreases the aw value of available substrate. Initially, in presence of substrate in the upper column, there is rise in total microbial count. However, with time, there is reduction of moisture which leads to the microbial growth in stationary phase, and with further reduction of moisture, the aw value reaches a point where the microbial growth cannot occur and the microbial growth curve enters the decline phase. Initially, there is rapid decrease in substrate by the high rate of consumption by growing microbes, but with time as moisture decreases continuously, microorganisms are unable to utilize it and the substrate does not decrease with lapse of time. Furthermore, the number of microbes reaching the higher depth is highly influenced by solid matrix of the soil column, which is indicated by decrease in total count of bacteria at higher depth. In addition, the amount of substrate appears to increase down the depth apparently due to reduction of moisture of soil.Figure 8
Variational trend of microbes in attached phase in depth and time.


## 5.1. Simulation Setup

Richards equation for unsaturated flow was linked to the microbial transport equation along with transport of growth-limiting substrate. In this simulation, we consider a vertical soil column of depthL=70 cm in a time period of tmax=1hr hr. The simulation starts with a uniform saturation θ=0.1 cm 3/cm3, and a constant water head ψ=−61.5 cm is maintained at the bottom boundary z=L. At the upper boundary z=0 (the soil surface), a constant flux qt=13.69 cm/hr for t<0.3 hr and zero normal flux condition for t>0.3 hr is adopted. Equal amounts of 103 cfu/ml microbial concentration in water phase and attached to solid matrix with 10 mg/100 ml of substrate concentration were set through column at 0.34 cm/sec of water flow velocity. Since no experimental evidence of the model is reported, all fixed model parameter values for soil characteristic used for model simulation are taken from Haverkamp et al. (1977). The values for the empirical parameters θs, θr, ks, ψ for soil type are adopted from [14]. In addition, the parameters related to the microbial transport and fate are listed in Table 1 ([6, 8]).Table 1
Kinetic parameters.
ParameterValueDescriptionD1.51×10−10m2/sDiffusion coefficientτ2.0/3.0Tortuosityvw0.34 cm/sSpeed of water flowv0.66 cm/sSpeed of microbe flowμm0.94 microbes/sSpecific growth rate of microbesks100 cm/sSaturation constantYs6.5 g/mol ATPBiomass produced per unit mass of substrate consumedB1.5×10−7microbes/hrDeath rate of bacteriaMma1100Maximum concentration of microorganisms in attached phaseka0.85Attachment coefficientkf0.95Freundlich isotherm slope parameterky0.37Detachment coefficientρs1Density of solid matrixThe moisture content has been taken as a function of the pressure head (suction); the constitutive relationship betweenθz,t and ψz,t is(45)θψ=θr+αθs−θrα+ψβ=0.075+1.611×106.287−0.0751.611×106+ψ3.96,Kψ=KsAA+ψγ=34×1.175×1061.175×106+ψ4.74.The simulation results are illustrated in Figures1–7. We show the variation trend of microbial concentration suspended in water Mz,t and attached to solid matrix Maz,t and the growth-limiting substrate Sz,t in depth.Figure 1
Front propagation analytic solution for the time evolution of the concentration of microbes in water phaseMz,t, plotted in relative units for several equidistantly separated depths z. The solution is obtained by solving (15) for M with θ given by the solution in (31), Ma approximated zero, γ=4.75, and the constants of material conveniently chosen. In the inset, we present the corresponding solutions σψ of (27) for several values for γ for comparison.Figure 2
Microbe concentration in waterMz,t (15) plotted vs. depth z for four moments of time. It can be seen from figure that the concentration of microbes renders due to the flow process and the consumption of growth limiting substrate but with decreasing in moisture content and because the property of aw microbial growth in stationary phase.Figure 3
Substrate concentration in waterSz,t plotted vs. depth z for four moments of time. It can be seen that the concentration of growth-limiting substrate decreases rapidly because of high consumption by microbes and increases apparently due to the decrease in moisture of the soil.Figure 4
Graphical representation of transport of microbes in aqueous phase for different depths (a) and substrates (b). Different colored curves represent distributions taken at moments of time separated by 0.1 hrs.
(a)(b)Figure 5
Graphical representation of transport of microbes in aqueous phase for different time periods (a), substrates, (b) and microbes in attached phase (c).
(a)(b)(c)Figure 6
Variational trend of microbes in aqueous phase in depth and time.Figure 7
Variational trend of substrate in aqueous phase in depth and time.The numerical procedure developed in the previous section is written in Python 3.7 Compiler with Intel® Core™ i5-2450M CPU @ 2.50 GHz.


## 5.2. Results and Discussion

The proposed model was applied to simulate the transport and growth of microbes with the transport of growth-limiting substrate through an unsaturated finite column. The growth of microorganisms is often controlled by a single substrate and the specific growth rate is assumed to be a function of the concentration of that rate limiting substrate. It also depends on the behavior of moisture, i.e., wetting and drainage, and the available substrate property, i.e., water activityaw of substrate. The aw of substrate is the ratio of vapor pressure of the substrate to the vapor pressure of pure water. Decrease in moisture content leads to decrease in aw value of the substrate that either decreases or ceases the microbial growth. The movement of microorganisms depends on the flow velocity, dispersion, diffusion, Brownian motion, and also chemostatic and tumbling motion of microbes. Here, in our mathematical model, we used the hypothetical values for microbes and growth-limiting substrate at the top of the soil column, and then the total microbial count, moisture content count, and substrate content count down the soil column were calculated. The model shows that the initial bacterial count, i.e., 1.1×103 cfu/ml, increases down the soil column in presence of substrate, i.e., 10 mg/100 ml, with the moisture of the porous medium. However, as the microbes and nutrients percolate down the column, there is steady decrease in moisture. The decrease in moisture showed direct influence on microbial growth, and the total microbial count decreased overall. Likewise, the substrate concentration showed decrease in concentration at the upper part of the column; however, the substrate appears to apparently increase towards the end due to reduction of moisture in soil. The decrease in microbial count down the column by continuous decrease in moisture results in adverse condition for microbial growth. Likewise, the apparent increase in substrate down the column is caused by decrease in moisture that concentrates the substrate in a limited amount of moisture. These phenomena are interpreted in the figures.Figure8 shows the front propagation analytic solution for the time evolution of the concentration of microbes in water phase. Figure 1 explains the growth of microbes in aqueous phase with respect to depth for flow velocity, whereas Figure 2 explains the growth of microorganisms attached to solid matrix with respect to depth for flow velocity. Furthermore, Figure 3 shows the consumption behavior of substrate. Figure 4 explains the transport of microbes and substrate for different time in aqueous and attached phase, and Figures 5–7 describe the variational trend of microbes and substrate in aqueous phase and attached phase in depth and time. The continuous decrease in moisture thus decreases the aw value of available substrate. Initially, in presence of substrate in the upper column, there is rise in total microbial count. However, with time, there is reduction of moisture which leads to the microbial growth in stationary phase, and with further reduction of moisture, the aw value reaches a point where the microbial growth cannot occur and the microbial growth curve enters the decline phase. Initially, there is rapid decrease in substrate by the high rate of consumption by growing microbes, but with time as moisture decreases continuously, microorganisms are unable to utilize it and the substrate does not decrease with lapse of time. Furthermore, the number of microbes reaching the higher depth is highly influenced by solid matrix of the soil column, which is indicated by decrease in total count of bacteria at higher depth. In addition, the amount of substrate appears to increase down the depth apparently due to reduction of moisture of soil.Figure 8
Variational trend of microbes in attached phase in depth and time.


## 6. Conclusion

In this work, we developed a mathematical model to predict the transport and growth of microorganisms in unsaturated porous medium (soil). The governing equation for microbial transport of microbes is coupled with another transport equation of substrate and the equation for the infiltration of water into a porous medium. Fully coupled numerical solutions of the governing equations have been obtained using finite difference and finite volume methods. Due to the stiffness of the Richards equation, we solved it using a stabilized explicit Runge–Kutta–Legendre super-time-stepping scheme [11]. For the transport equations, being hyperbolic in nature, we employed an upwind finite volume method for the advection terms. The numerical simulations offer realistic solutions to the mathematical model as shown in Figures 1–4 and thus open up the possibilities of prediction and understanding of the described phenomena and comparison with the experimental data. Following the auxiliary equation method for nonlinear diffusion-reaction equations, we also obtained a closed form traveling wave solution for the Richards equation. The closed form analytic solution contains several parameters (constants of integration). Variation in these parameters corresponds to variation in the physical constants (material property) present in the mathematical model, thereby representing change in different experimental setup. After many hit-and-trial experiments, we obtained a traveling wave solution with a specific wave front velocity which mimics the qualitative behavior of the solution of the nonlinear system. The authors believe that the occurrence of these (analytically obtained and numerically validated) front moving solitons for the transport model raises an interesting problem for the researchers to look at the solution using the data science/machine learning approach.



---

*Source: 6278126-2021-12-27.xml*