Numerical modeling of virus transport through unsaturated porous media

: This paper describes the movement of virus in one-dimensional unsaturated porous media. The governing virus transport equations consider the inactivation in liquid phase, liquid–solid interface, air–liquid interface, and sorption in both liquid–solid and air–liquid interfaces. Finite-volume method has been used for solving the advection and dispersion processes of the virus transport equation. The effects of transport parameters on virus concentration profiles have been investigated for virus present in liquid phase, adsorbed liquid–solid and liquid–air phases. The results show that the movement of viruses in three phases is affected by soil moisture, inactivation rate, pore velocity, and mass transfer coefficients. It is found that the magnitude of virus sorption is higher at the air–liquid interface as compared to the liquid–solid interface. A higher value of mass transfer coefficient leads to an increase in the virus concentration in both liquid–solid and air–liquid interfaces.


Introduction
In the past few decades, it has been noticed that the level of groundwater contamination has increased due to increase in industrial and agricultural activities. The presence of viruses in drinking water causes human diseases and it originates from septic tanks, sewage sludge, sanitary landfills, and agricultural practices. The experimental studies reported in the past indicate that the virus survives for a certain period of time in unsaturated porous media before reaching into the subsurface water (Chu, Jin, Baumann, & Yates, 2003;Schaub & Sorber, 1977;Yates & Ouyang, 1992;Yates, Yates, Wagner, & Gerba, 1987). Hence, it is essential to understand the transport process of viruses through the unsaturated porous media to prevent further contamination of subsurface aquifer system.

PUBLIC INTEREST STATEMENT
In this study, finite-volume method has been used to get the numerical solution of virus transport equations considering inactivation in liquid phase, liquid-solid interface, air-liquid interface, and sorption in both liquid-solid and air-liquid interfaces. The present model has been used to investigate the behavior of virus concentration in liquid phase, adsorbed in liquid-solid and air-liquid interfaces with depth. Present model can be used to simulate the experimental data of viruses in the field and also to predict the movement of viruses in the subsurface media.
Unsaturated porous media consist of liquid, solid, and air phases and for water-wet solid surfaces, both liquid-solid and air-liquid interfaces exist (Freeze & Cherry, 1979). It is also known that the unsaturated zone plays an important role in the transport of fluids and contamination from the surface to the groundwater. Virus sorption within unsaturated porous media is affected by the presence of two interfaces. Viruses are sorbed onto liquid-solid interfaces via physical adsorption, chemical sorption, or ion exchange (Vilker & Burge, 1980). Vilker (1981) has suggested that non-equilibrium kinetic sorption is appropriate for describing virus attachment to the liquid-solid interfaces and for viruses with size similar to the size of solutes. This sorption process represents the rate of approach to equilibrium between adsorbed and liquid-phase virus concentration considering virus transport to the outer layer of a solid particle by mass transfer. Powelson, Simpson, and Gerba (1990) have suggested that the virus sorption is greater at air-liquid than liquid-solid interfaces. Virus sorption takes place at the liquid-solid interfaces due to electrostatic double-layer interactions and van der Waals forces (Teutsch, Herbold-Paschke, Tougianidou, Hahn, & Botzenhart, 1991). Tim and Mostaghimi (1991) developed a numerical model for water flow and virus transport in variably saturated porous media assuming that the virus sorption is an equilibrium process. Park, Blanford, and Huyakorn (1992) developed a semi-analytical model for both steady state and transient vertical virus transport in the unsaturated media and along the flow lines in the saturated zone considering equilibrium sorption and first-order inactivation. Yates and Ouyang (1992) developed a one-dimensional numerical model for flow of water, viruses, and heat in unsaturated porous media considering moisture-independent sorption, filtration, and temperature-dependent inactivation. Adsorption and inactivation are two different processes of virus removal, and viruses can get detached from the soil media because of the ions present in the groundwater (Bales, Li, Maguire, Yahya, & Gerba, 1993). Virus sorption at air-liquid interfaces is controlled by virus particle surface hydrophobicity, solution ionic strength, and particle charge (Wan & Wilson, 1994). Poletika, Jury, and Yates (1995) have showed that the viruses at an air-liquid interface may be desorbed under high interfacial shear stresses induced by fast interstitial fluid flow.
Virus inactivation is generally considered as a first-order irreversible sink mechanism (Sim & Chrysikopoulos, 1996). Recently, the mathematical models have been developed for virus transport in saturated porous media considering different inactivation rates for viruses in different phases (Sim & Chrysikopoulos, 1996, 1999. Schijven and Šimůnek (2002) have shown that the factors like size, attachment characteristics, and rate of inactivation affect the movement of viruses in porous media. Due to advection and dispersion, viruses get spread in the soil media which reduces its concentrations. The removal of viruses from the groundwater occurs due to the processes of adsorption and inactivation (Chattopadhyay, Chattopadhyay, Lyon, & Wilson, 2002;Chu, Jin, Flury, & Yates, 2001). Various factors like temperature, moisture content, pH, hydraulic conditions affect the transport of viruses below the ground surface. However, the temperature significantly affects the transport of viruses and a relationship between increase in temperature and inactivation rate of virus for different viruses has been given by Gerba and Rose (2003), and Gerba and Smith (2005). Torkzaban, Hassanizadeh, Schijven, de Bruin, and de Roda Husman (2006) studied the transport of bacteria through saturated and unsaturated porous media and their results demonstrate that the attachment to the air-water interfaces is reversible. Anders and Chrysikopoulos (2009) conducted soil column experiments under both saturated and unsaturated conditions. Their results indicate that even for unfavorable attachment conditions within a sand column, saturation levels can affect the virus transport through porous media. Further, many researchers have used the numerical method to investigate movement of viruses through saturated and unsaturated porous media and also estimated virus transport parameters considering liquid and solid phases (Joshi, Ojha, Sharma, & Surampalli, 2013;Ratha, Prasad, & Ojha, 2009;Sharma & Srivastava, 2011). Syngouna and Chrysikopoulos (2015) studied the effect of colloids and water saturation level on the attenuation and transport of colloids and viruses in unsaturated porous media. Thus, from the past research works, it is found that the effects of transport parameters on viruses in three phases, i.e. liquid phase, adsorbed liquid-solid and liquid-air interfaces have not been studied in detail. Hence, the present study describes the finite-volume method approach to solve the governing transport equation and investigates the behavior of movement of viruses in three phases through unsaturated porous media.
In this study, the numerical finite-volume method has been used to solve the governing equations for virus transport in one-dimensional, unsaturated porous media considering virus sorption on liquid-solid and air-liquid interfaces. An attempt has been made to investigate the effects of moisture content, inactivation rate constant, pore velocity, and mass transfer coefficients on the variation of virus transport along with vertical depth.

Virus transport
The governing differential equation for virus transport in one-dimensional, homogeneous, unsaturated porous media can be written as (Sim & Chrysikopoulos, 2000): where C represents the virus concentrations in liquid phase (ML −3 ); C s represents the adsorbed at the liquid-solid interface (MM −1 ); C a represents the adsorbed at the air-liquid interface (ML −3 ); θ represents the volumetric moisture content (L 3 L −3 ); q represents the specific discharge (LT −1 ); , s , and a are the inactivation rate constants of viruses in the liquid phase, adsorbed in the liquid-solid interface and adsorbed in the air-liquid interfaces (T −1 ), respectively; and ρ is the bulk density of the soil media (ML −3 ). The hydrodynamic dispersion coefficient D z can be expressed as (L 2 T −1 ) (Nielsen, Van Genuchten, & Biggar, 1986): where α z is the dispersivity (L), and D 0 is the molecular diffusion coefficient (L 2 T −1 ).

Virus sorption at interfaces
The expressions for the viruses adsorbed at solid-liquid and air-liquid interfaces can be written as (Sim & Chrysikopoulos, 1999): where k represents the mass transfer rate of liquid to liquid-solid (T −1 ); k a liquid to air-liquid interface (T −1 ); and C g is the virus concentration in liquid phase (ML −3 ), which is in close contact with soil solids. It was assumed a linear equilibrium relationship by Sim and Chrysikopoulos (1996) and it can be expressed as: where K d is the distribution coefficient (L 3 M −1 ); κ is the mass transfer coefficient for liquid to liquidsolid interface (LT −1 ); and a T is the specific area of liquid-solid interface (L 2 L −3 ). It can be defined as the ratio of total surface area of soil particles to the bulk volume of the porous medium and it is expressed as (Fogler, 1992): where r p is the average radius of soil particles (L); and θ s is the water content of a saturated porous medium. The mass transfer rate coefficient (k a ) for liquid to air-liquid interface (T −1 ) is given as: where κ a is the mass transfer coefficient of liquid to air-liquid interface (LT −1 ); a Ta is the specific area of air-liquid interface (L 2 L −3 ); and the expression can be given as (Cary, 1994): The expression for r 0 , i.e. the effective pore radius can be expressed as: where ζ and b are the soil constants; σ is the surface tension of water (MT −2 ); ρ w is the density of water (ML −3 ); g is the gravitational constant (LT −2 );and h 0 is the air-entry value (L), which is defined as the pore water head where air begins to enter water-saturated pores. It can be observed from Equation (9) that a Ta is a function of moisture content and its value decreases as the available moisture content increases and becomes zero when the available moisture content is equal to the saturated moisture content, i.e. θ = θ sat .

Initial and boundary conditions for virus transport
Initially at time t = 0, it is assumed that there is a negligible concentration of viruses present in all the three phases of soil media.
where C 0 is the source concentration of virus on the ground surface.

Numerical model
Numerical method admits any arbitrary boundary condition and it can be used in cases of dealing with complex problem for which analytical solution cannot be obtained. The numerical method can also be used for real-field problem such as complex geometry, partial variation of hydraulic conductivity and non-linear problem. However, analytical method has limitations in case of non-linear problem; in such cases, a numerical method works very well.
The finite-volume method has been used to get the solution of governing equations for virus transport. Using Equations (3) and (4), the simplified form of Equation (1) can be expressed as: The advective term of Equation (12) can be given as (Putti, Yeh, & Mulder, 1990): The dispersion transport term can be expressed as: The operator split approach reduces the reaction terms in Equation (12) to a coupled ordinary differential equation which can be expressed as:

Procedure of numerical solution
The finite-volume method is used for solving the advective transport and is based on monotone upwind schemes for conservation laws (MUSCL) (Van Leer, 1977). Equation (12) has been broken into three parts, based on the suggestion by Putti et al. (1990), in which the explicit numerical scheme is used for advective transport and an implicit numerical scheme is used for dispersive transport, and explicit method for reaction part. The advantage of this method is to handle either advection dominated or dispersion dominated for solute transport through porous media accurately. This method is globally high-order accurate and non-oscillatory, and the detailed procedure for solution of the advective part has been given by Ratha et al. (2009), and Sharma, Joshi, Srivastava, and Ojha (2014).

Dispersive transport equation
The resulting output concentrations of the advective transport are used as the initial condition for dispersive transport. A conventional fully implicit finite-difference scheme, which is unconditionally stable, is used to obtain the final concentration at the end of time step. The implicit finite-difference formulation of dispersive transport equation can be expressed as: where Δz and Δt are spatial grid size and time step, respectively.
Remaining Equations (3), (4), and (15) are solved using the explicit numerical method and the formulation is given below:

Verification of model
The present numerical model is verified using the analytical solution given by Sim and Chrysikopoulos (1995). The numerical and analytical results of spatial virus concentration in liquid phase have been shown in Figure 1 with two different values of inactivation rate constants. The grid size Δz = 0.2 cm and time step Δt = 0.025 h have been used in this simulation. It can be seen that the numerical results match very well with the analytical solution. To reduce the numerical error, the values of both the Peclet number and Courant number are kept less than one.

Results and discussion
The numerical model has been used to investigate the behavior of virus concentration profiles through unsaturated porous media. The input parameters used for simulation are shown in Table 1 (Sim & Chrysikopoulos, 2000), and a continuous source of virus is injected at the ground surface. As the moisture content of porous media reduces from its full saturation to partial saturation, there is a significant decrease in the liquid-phase virus concentration and significant increase in the virus concentration sorbed at air-liquid interface. This increase is due to the increase in specific air-liquid interface.
Figures 3a-3c show the variation of virus concentration along with depth for different values of inactivation rate constants. In this study, the inactivation rate coefficients λ, λ a , and λ s all are assumed to be constant. However, in reality, its value depends on temperature and time (Sim & Chrysikopoulos, 1996). It is assumed that λ a = λ and λ s = λ/2 as suggested by Yates and Ouyang (1992), and Thompson, Flury, Yates, and Jury (1998). The simulation is carried out with constant pore velocity v = 4 cm h −1 , uniform water content θ = 0.35, and saturated water content θ sat = 0.45. As expected, with an increase in the value of inactivation rate coefficients, the concentration of viruses in liquid phase, liquid-solid interface and air-liquid interface decreases. The behavior of virus concentration is non-uniform along the depth.    -5c show the effect of mass transfer coefficient on the variation of concentration of viruses in the liquid phase, liquid-solid interface and air-liquid interface. The simulations are carried out with constant pore velocity of 4 cm h −1 , considering negligible inactivation constants. It is assumed that the mass transfer coefficient of liquid to air-liquid interface is five times the mass transfer coefficient of liquid to liquid-solid interface as considered by Sim and Chrysikopoulos (2000). Three different values of liquid to liquid-solid mass transfer coefficients, i.e. κ = 0.006, 0.012, and 0.024 cm h −1 . It is seen that on increasing the value of mass transfer coefficient (κ), the concentration of viruses in liquid phase decreases and concentration of viruses in both liquid-solid interface and air-liquid interface increases. The magnitude of virus concentration in air-liquid interface (C a ) is higher in comparison to the concentration of viruses in liquid-solid interface (C s ) for the same depth. This has occurred due to the assumption of k a = 5 k.

Summary and conclusions
In this study, a finite-volume method has been used to develop a numerical model to analyze the virus transport in one-dimensional unsaturated porous media. The model accounts for the virus sorption on the liquid-solid and air-liquid interfaces as well as inactivation of viruses suspended in the liquid phase and virus attached to both interfaces. The effects of transport parameters on virus concentration profiles were investigated and the findings are listed below: (1) In the field, the moisture content of porous media reduces from its full saturation to partial saturation; there is a significant decrease in concentration of viruses in liquid phase and liquid-solid interface, and an increase in concentration of viruses in air-liquid interface. This increase is due to the increase in liquid to air-liquid mass transfer rate, as it depends on the specific air-liquid interface area, which increases with a decrease in the value of moisture content.
(2) As expected, with an increase in the value of inactivation rate coefficients, the magnitude of virus concentration in liquid phase and both adsorbed interfaces reduce. A higher value of pore water velocity leads to increase in the movement of viruses in depth. It is also seen that the variation of virus concentration along the depth is non-uniform.
(3) A higher value of mass transfer coefficient decreases the magnitude of virus concentration in liquid phase, while it increases virus concentration in both liquid-solid and air-liquid interfaces. Also, the magnitude of viruses in air-liquid interface is higher in comparison to concentration of virus in liquid-solid interface at the same depth. Finally, this model may be used for simulation of experimental data of viruses in the field.