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.
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 Section 2, 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 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 . If one assumes gravity () as the only driving force of the flow, one maps 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 Section 3, 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 and the same thing for the pressure head 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 direction ( downward) with no source/sink term has the formwhere is the volumetric moisture content, is the pressure head, and 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]:where and represent the saturated and residual moisture content, respectively; corresponds to the saturated hydraulic conductivity; and , , , 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].
Let denote the concentration of microbes in the water phase defined as the mass of microbial cells per unit volume of water, the bulk density of the porous medium, and the mass of microbes attached to the solid phase. At moment , the total biomass in the elementary volume is the sum of the mass in the water phase plus the mass in the attached phase [6, 8]. The governing equation for microbial transport in the porous media in one dimension is given by the following partial differential equation:where is the velocity of water flow, is the chemostatic velocity of the microbes, is the tortuosity of the pore space, is a diffusion coefficient, is the growth rate of microbes combined for water and attached biomass phase (mass of microbes per unit volume and time), and 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 aswhere is the mass of substrate per unit volume of water, is the mass of substrate adsorbed per unit mass of solid matrix, is the diffusion-dispersion coefficient of the substrate, and 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 , , , and 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 aswhere is the specific growth rate of microbes given in [15] in the formwhere is specific growth rate of microbes and is saturation constant. The decay rate of microbes is defined aswhere is the decay or inactivation constant. The substrate consumption is given bywhere 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 variables in (2.4) and in (4). For the substrate in the attached biomass, we assume equilibrium adsorption isotherm (Freundlich) [16] aswhere is Freundlich isotherm slope parameter. For the microbial growth in the attached biomass phase, we use the following nonequilibrium microbial kinetics:where is the rate of microbial attachment and is the rate of microbial detachment. Following [6], we definewhere is the maximum concentration of microbes in attached phase, is attachment coefficient, and is detachment coefficient.
2.5. Mathematical Model
By taking 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 nonlinear PDE system in the functions , and in a bounded region as follows:subject to the following initial and boundary conditions:(i)Initial condition: , , , .(ii)Lower boundary condition: , , .(iii)Upper boundary condition: 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 . 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 formwhere we consider and to be invertible differentiable functions, at least for conveniently chosen intervals for the variables, and where we denote
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 and , (17) becomes a Burgers equation:in , 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 conditionwhere we assume the function is invertible, differentiable, and nonzero and where is an arbitrary constant. With this type of mapping , (17) becomeswhich is a diffusion equation with nonlinearity resulting from its functional coefficients. For example, if is given by the material relation in (3), by integration of (20), we havewhich is a strictly monotonic function. For small , (22) generates a linear dependence , and for larger , (22) approaches a 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 approximationwe obtain for (21) an exact traveling front solution [18]:where is the comoving coordinate; is a free parameter representing the velocity of the front; and are constants of integration. For example, if we choose a simple case where , the following results:
Consequently, in the first order of their Taylor series, we have
In continuation, we can choose the arbitrary parameters 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 which would represent the choice of a certain value for the pressure head generating maximum, or minimum, moisture content. Consequently, which means a choice . In this case, the solution equation (24) of (21) has an even simple form:where is constant of integration and Ei is the exponential integral function. Because the inverse of the exponential integral function is multivalued, the solutions 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 , (17) can be written in the comoving frame in the formwith
We can write the differential equation fulfilled by the inverse function in the formwhich has solutionwhereand the traveling solution 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 . We consider in a first approximation that , and we take a realistic value for (see Section 5) and for the other constants of material. In the comoving frame , (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:where the functional coefficient is given bywhich depends on the pressure head through the Kirchhoff transformation defined aswith and . 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 length ) into finite cells , with length and define the nodes as the midpoint of the cell . That is, the nodes and the faces of the control volumes are given by
Let be the time increments and define the discrete time-steps , .
The initial time is denoted by , and the discrete approximation for at the grid is denoted by . We also regard as an approximation to mean value of in the cell at time as
The variables and can be derived in similar manner (13)
If an approximate solution is assumed to be known at all grid points at time , a method must be specified to advance the solution to time , subject to the boundary conditions. For this, in finite volume method, we integrate equation (14) into the finite cell and approximate the integral to get the numerical scheme.
The first term in the integral (38) is approximated using forward difference in time as
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 . Clearly, these fluxes do have two parts, advection and diffusion. For the advective fluxes we use upwind approximation.
Approximating the third term asyields the following scheme:
Similarly, we write the finite volume discretization for (14) and (15).
Discretization of the growth model equation (14) is as follows:
Discretization of the microbial transport equation (15) is as follows:
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 depth cm in a time period of hr. The simulation starts with a uniform saturation cm , and a constant water head cm is maintained at the bottom boundary . At the upper boundary (the soil surface), a constant flux cm/hr for hr and zero normal flux condition for hr is adopted. Equal amounts of 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 , , , 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]).
The moisture content has been taken as a function of the pressure head (suction); the constitutive relationship between and is
The simulation results are illustrated in Figures 1–7. We show the variation trend of microbial concentration suspended in water and attached to solid matrix and the growth-limiting substrate in depth.
(a)
(b)
(a)
(b)
(c)
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 activity of substrate. The 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 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., 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.
Figure 8 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 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 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.
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.
Data Availability
The data used for supporting the findings of this study are included within the article.
Conflicts of Interest
The authors declare that they have no conflicts of interest.