Simulation of Contaminant Transport through the Vadose Zone: A Continuum Mechanical Approach within the Framework of the Extended Theory of Porous Media (eTPM)

: The simulation of contaminant transport through the vadose zone enjoys high signiﬁcance for decision makers and contaminated site planners since the vadose zone can serve as a ﬁlter, but many contaminants can be transported from this region to aquifers. The intention of this paper is to utilize the extended Theory of Porous Media (eTPM) to develop a ternary model for the simulation of contaminant transport in the vadose zone whose application is subsequently shown via a numerical example. The simulation was conducted for 140 days, during which the contamination source was removed after 25 days. The results indicate that the contaminant reached the water table after 76 days. The concentration of the contaminant reaching the groundwater was 17% less than that of the contaminant source.


Introduction
The vadose zone is the connection point between the surface of the earth, precipitation, surface water, and aquifers and enjoys great significance in the hydrological cycle, including water and energy in the land surface system of the earth [1,2].This zone adjusts the water flow from the land surface to the groundwater and vice versa, and it plays a crucial role in the transport of solutes and their retention.The most severe contamination threats to groundwater quality emanate from coincidental surface releases of various types of wastes and products associated with industrial, agricultural, or urban development, such as leachate from municipal landfills [3].These contaminants are transferred through the vadose zone, which is partially saturated with air and water.Although the vadose zone functions as a buffer and filter zone, in which attenuation, absorption, screening, degradation as well as chemical and biological transformations of contaminants occur, most of the hazardous contaminants are reactive in soil, and they infiltrate through the vadose zone and permeate unsaturated soil in order to reach the aquifer and continue to migrate in the direction of the groundwater flow [4].Hence, the contamination of the vadose zone can be considered as a continuous source of contamination for aquifers, and it affects the quality and quantity of groundwater resources [5,6], leading to significant risks for human health and the environment.
The risk assessment, which is required by decision makers for the long-term management of the contaminated vadose zone, necessitates an understanding of the water flow and the fate of contaminants and their concentration at certain points [7].Because the direct measurement and monitoring of contaminants' concentration and the impact of a contaminated vadose zone on the groundwater quality are difficult and costly processes in many cases, mathematical modelling methods can instead be employed to quantify such impacts.Furthermore, the results of such models can be utilised for testing alternative conceptual models, predicting future conditions, providing information required for remediation design, interpreting and analysing field experimental data, and subsequently evaluating corrective measures.
In view of this, several mathematical models have been developed based on various analytical and numerical studies.In general, there are two different approaches to the simulation of flow and contaminant transport in the vadose zone, namely analytical [8][9][10] and numerical methods [11,12].Due to the complexity of phase partitioning and reaction, there are only a few analytical models considering multiple phases and multiple solutes [13,14], multiple solute transport with reactions in multiple phases [15], and combined transient water and solute transport [16].In order to deal with problems concerning saturated and unsaturated zones, there are different numerical approaches, including the finite difference method (FDM) [17,18], the finite volume method (FVM) [19], the finite element method (FEM) [20,21], and the Meshfree method [22], each of which has various treats that tailor the model to the special application for which it was designed.Rao et al. combined explicit finite difference formulations with multiple domain algorithms and developed models for one and two-dimensional transient contaminant transport flow [23,24].Kacur et al. introduced a continuum mechanical model for fluid flow in unsaturated deformable porous media, where the fluid could undergo phase transitions [12].Kuntz et al. developed a numerical simulation of the reactive transport of compounds in the vadose zone to investigate whether steady-state flow simulations accurately describe reactive contaminant transport under transient conditions in the field [25].Bunsri et al. developed a model for water movement and contaminant transport through unsaturated soil based on Richard's equation and utilizing the finite element method [26].Rajsekhar et al. described virus transport in partially saturated soil using the finite volume method; they considered the virus sorption on the liquid-solid and air-liquid interfaces as well as the inactivation of viruses suspended in the liquid phase and viruses attached to both interfaces [27].Furthermore, different software packages, including HYDRUS, Modeflow, and Chemflux, have been developed over the last few years, each having its own advantages and disadvantages [28,29].HYDRUS describes the vadose zone process using the one-dimensional Richards equation, which requires moderate computational effort [30].
Most of the flow and contaminant transport models of the vadose zone do not consider the entire coupling of the thermo-hydro-mechanical interaction, phase transitions, deformability of the soil, physical phenomena concurrent with water flow, air flow, heat transfer, and contaminant transport, and in particular, the effects of air flow, temperature, and chemical reactions.The Theory of Porous Media (TPM) [31][32][33][34] facilitates a comprehensive continuum mechanical framework based on the Mixture Theory [35], enhanced by the concept of volume fractions, for the simulation of different processes in a saturated or partially saturated medium.Furthermore, it paves the way for the description of geophysical and geochemical environmental problems [18,21,[36][37][38][39][40][41][42].
This study aims to provide a numerical simulation of contaminant transport in the vadose zone within the framework of the extended Theory of Porous Media.In developing the simulation, the vadose zone is modeled as triphasic system consisting of a soil phase (S), saturated by water as a liquid phase (L) and a gaseous air phase (G).The phases are assumed to be immiscible and statistically and homogeneously distributed over the control volume.Each phase can be composed of or can contain miscible components, which can be considered as contaminants.Their motion through porous media in regard to different transport regimes, such as diffusion and advection, can be simulated.
In this paper, the development of a triphasic model based on the Theory of Porous Media (TPM) is initially presented.The model is expanded by considering a contaminant as a solute in the water utilizing the extended Theory of Porous Media (eTPM).For the conceptual and numerical testing of the model, a numerical example is considered and its simulation results are presented.

Theory of Porous Media (TPM)
The fundamentals of the TPM were developed from the combination of two basic concepts, namely the Theory of Mixtures [43] and the concept of volume fractions [44][45][46], which enable the consideration of the local composition of the overall aggregate.The vadose zone as a partially saturated soil, denoted as ϕ, is composed of phases ϕ α with α ∈ {S, L, G}, namely pure soil as a solid (S) matrix, which generates the control volume B S and its boundary ∂B S , saturated by a fluid (F) which consists of liquid (L) and gas (G), expressed as follows: ( The TPM utilizing the Theory of Mixtures smears each heterogeneously distributed constituent of the true structure over a control volume.The applied homogenization procedure of each phase of the vadose zone over a representative elementary volume (REV) leads to the ternary superimposed continuum model, see Figure 1.As a consequence of the superposition of the homogenized phases, each phase ϕ α occupies each spatial point of the control volume simultaneously.The volume fraction n α of each phase is defined with where x denotes the position vector in the current configuration at time t, dv α is the partial element volume of the phase ϕ α , and dv is the element volume of the bulk medium [32,38].
With the definition of the volume fraction, which assigns a volumetric dimension to each phase in every spatial point x, each phase can be identified within the mixture.Due to the applied homogenization, the partial density ρ α is defined as the partial mass dm α of phase ϕ α divided by the bulk volume dv, expressed by With Equation ( 2), the partial and realistic density ρ αR are connected via Regarding the constraint that no vacant space within the control volume B S is allowed, and with the result that the solid matrix is entirely saturated by the fluid phases water and air, the saturation condition can be written as where x denotes the position vector in the actual configuration of the material point X α .Balance equations describe the equilibrium condition that each natural system strives for.Applied mechanics deals with four balance equations, which draw a conclusion from the mass, the momentum, the moment of momentum, and the energy of a continuous body.Analogously to a single-component material, these balance equations can be written for all phases in the framework of the TPM, considering all interactions for the description of chemical and physical exchange processes, such as the mass or energy exchange between the phases.The balance equations in the aforementioned order can be written as follows: The quantities ρα , pα , êα , ε α , r α , q α represent mass, momentum, and energy supply, as well as the specific internal energy, an internal heat source, and the heat flux over the surface, respectively.Moreover, b, D α and T α denote the body force vector, the symmetric part of the spatial velocity gradient, and the partial Cauchy stress tensor, respectively.The detailed kinematic relations are described in Appendix A. Since the sum of the local balance equations over all phases must comply with balance equations for a one-component material and for the whole mixture, the sum of the supply terms has to be zero [47]: The specific internal energy ε α can be rewritten in terms of the HELMHOLTZ free energy function ψ α and the specific entropy η α with the following equation: where θ α is the absolute phase temperature.Establishing the common entropy inequality for all phases is a necessary and sufficient condition for the existence of a dissipation mechanism within the mixture [32] and provides thermodynamically consistent restrictions for constitutive modelling.Hence, the entropy inequality must be evaluated for the whole mixture by the sum over all κ phases with the following equation: where 'div' denotes the spatial divergence operator.

Extended Theory of Porous Media (eTPM)
To represent miscible substances within the macroscopic immiscible phases ϕ α , the TPM has to be extended [48].The extended TPM (eTPM) offers a comprehensive approach for the simulation of different physical and chemical processes, such as the transport processes of miscible substances, i.e., solutes, in their carrier phases (solvents).On the basis of the eTPM, the composition of the porous medium, which is the object of our investigation, can be rewritten as where ϕ αβ denotes the ν miscible substances ϕ β solved in the carrier phase ϕ α .To consider the contribution of the solute to the whole mixture, the mixture molar concentration and the solvent molar concentration of the solute can be defined as follows: c where dn β mol denotes the number of moles, while dm β and M β mol define the local mass of a solute ϕ β in terms of the molar mass.Analogously to Equation (3) and Equation ( 4), different mass densities can be defined for the miscible components.The mass concentration of component ρ β in the solvent ρ α is given with and its partial density as The mass and molar concentration (Equations ( 16) and ( 15)) are connected via the molar mass with In analogy to the macroscopic and immiscible quantities, the balance equations, kinematics, and constraints of the supply terms are valid for the mixture components.

Model Assumptions
The vadose zone is composed of the soil matrix ϕ S saturated by water, represented by the liquid phase ϕ L , and by air, represented by the gas phase ϕ G .The liquid phase is composed of the contaminant ϕ LC solved in the carrier phase ϕ L , so that the mixture body under investigation can be described with Since the volume fraction of the water contaminant is negligible in comparison to that of the carrier phases, no volume change dependent on a concentration variation is considered, and the saturation condition is given as follows: For the description of the biphasic fluid mixture, it is convenient to introduce the saturation function s β , which defines the ratio of the pore fluids in the available pore space with where n F = n L + n G denotes the porosity variable, which leads to Therein, s L and s G denote the saturation degrees of the liquid and gas phases, respectively.
Soil and water are considered to be incompressible; thus, their realistic densities ρ SR and ρ LR are constant, resulting in ρ SR S = ρ LR L = 0.The air is assumed to behave as an ideal gas, and its density change depends on the temperature θ and the effective air pressure p LR , with ρ LR = ρ LR θ, p LR .Mass transfer is not taken into consideration, such that The constraint for the interaction forces (Equation ( 10) 2 ) evaluated for the investigation of the contaminated vadose zone reads with All constituents investigated contain the same temperature with θ S = θ L = θ L .Thus, there is no energy supply êα between the phases, such that Equation (18) 3 vanishes with êS = êL = êG = 0. ( Additionally, both a constant temperature and no internal heat supply source r α are considered, so that the energy balance within the set of governing equations can be excluded.Moreover, quasi-static conditions are assumed for all macroscopic phases as well as for the contaminant, so that the acceleration terms can be neglected.Although this restricts the validity of the model to slow processes, it represents a reasonable assumption since acceleration terms fundamentally enact no part for the vadose zone, except in some circumstances such as during an earthquake.In addition, the square of the velocities x α • x α and x LC • x LC are not taken into consideration.

Field Equations
Considering the introduced assumptions, the quasi-static, triphasic, multi-component model is derived to describe water flow and contaminant transport in the vadose zone.A set of coupled field equations, balance equations of mass and momentum, is required to calculate the unknown quantities while considering the restrictions associated with the volume fractions and interaction forces.The local forms of the mass balance equations for the solid, liquid, and gas phase are given with where w LS = x L − x S and w GS = x G − x S denote the relative difference velocities between the liquid and gas phase, and the solid deformation.The local mass balance equation of the contaminant is derived as Considering the constant molar mass M C mol and the time derivative with respect to the solid phase leads to where j LC mol = n F s L c LC mol w LCS indicates the molar flux of the contaminant, which is given with the difference velocity w LCS = x LC − x S .Considering the given assumptions, the partial local momentum balance equations for the individual phases and the contaminant read as follows:

Constitutive Model
In addition to the kinematics and balance equations, constitutive equations are required in order to complete the system of equations and to identify and define the partial Cauchy stress tensors T S , T L , T G and T LC , the linear momentum productions pG , pL , and pLC , the saturation s L and the effective gas density ρ GR .Proceeding from general thermodynamical considerations [49,50] and the evaluation of the entropy inequality (see Appendix B), the partial Cauchy stresses are given with the following relations: In the equation above, K S = 1 2 (B S − I) denotes the Karni-Reiner strain tensor, with B S = F S F T S being the left Cauchy-Green tensor.The quantities p LR , p GR , and p FR denote the effective pressure of the wetting pore water, the effective pressure of the non-wetting pore air, and the effective overall fluid pressure, respectively, and µ LC denotes the chemical potential of the contaminant.
Incorporating the momentum production terms into the corresponding momentum balances, the individual relation for the filter velocities, which can be defined as the product of the volume fraction (see Appendix A.2), are given as follows: where µ L and µ L indicate the shear viscosity of air and water, K S is intrinsic soil permeability, and k G r and k L r are the relative permeability of air and water.These relations represent the standard Darcy law derived from the experiments for the multi-phase flow problems [51].
In addition, the flux of the contaminant in the water can be established with where D Fick , R and θ denote the diffusivity of the contaminant within the water, universal gas constant and absolute temperature.

Numerical Treatment
By utilizing the FEM in the framework of a standard Galerkin procedure, the weak formulation of the local balance equations is formulated.The transformation of the strong local form of the mechanical balance equations into a weak global representation is realised by multiplying them with the following independent weighting functions: In order to derive the weak formulation of the momentum balance of the overall mixture, this equation is multiplied scalar by the respective test function δu S representing a virtual deformation of the solid phase.Integrating over the domain B S and restating the term governed by the divergence operator by utilizing the product rule and the Gaußian integral theorem yields the weak formulation as Therein, T LC is neglected (cf.Appendix A.1).The surface integral appears, which is governed by the overall surface traction t = T S + T L + T G n acting on the so-called Neumann boundary ∂B S of the overall mixture.Analogously to the overall momentum balance, the weak forms of the volume balances of the water and air can be derived by weighting them with the test functions δp LR and δp GR , respectively.This yields and where n L w LS • n and n G ρ GR w GS • n denote the flux of the water volume and air mass through the Neumann boundary ∂B S .The weak form of the concentration balance for the solute is derived with

Results and Discussion
The simulation of contaminant fate in the vadose zone is an ongoing research topic, wherein a variety of numerical studies are performed to determine the dynamics of the vadose zone.The computations presented here resulted from a triphasic TPM formulation; they were used together with the numerical tools presented in the previous sections and implemented into the finite element program FEAP, developed by Taylor, UC Berkeley.To ensure stable numerical results, the so-called Taylor-Hood elements were used with quadratic ansatz functions for the displacements, while the linear ansatz functions were used for the remaining degrees of freedom.The triphasic model was applied to a numerical example.Table 1 illustrates the material parameters and initial values, which were used to describe the vadose zone and contaminant transport.
The numerical example corresponds to a two-dimensional description of the contaminant transport through the vadose zone (cf. Figure 2).The vadose zone has a height of 40 m and a width of 5 m.Furthermore, the left side of the vadose zone is assumed to be more or less impermeable.In addition, the right side of the vadose zone is impermeable; only the last 10 m on the right side are permeable, and the filter on the right-hand side of the structure is assumed to be governed by K S = 10 −8 m.The vadose zone is governed by an intrinsic permeability of K S = 10 −10 m.Based on the triphasic formulation, it is assumed that a water table on the left side of the vadose zone is slowly increased.As a result, water streams into the vadose zone and finally reaches a stationary state.

Numerical Example
The level of the groundwater table can vary because of natural phenomena such as precipitation and oceanic oscillations such as tides and waves in the coastal aquifers [52]; the variation can also be caused by human activities.The change in the groundwater table is determinative of the variation of the saturation degree, seepage velocity, and contaminant fate.In order to develop a remediation process, it can be beneficial to know how varying the groundwater table and saturation degree in the vadose zone can affect the contaminant fate.Concerning this, the sketched vadose zone (see Figure 2) is numerically investigated together with the changing of the groundwater table on its left side.The example shows a two-dimensional description of the contaminant transport through the vadose zone to the water table.The investigated vadose zone has a dimension of 5.0 m x 40.0 m.The simulation was conducted for 140 days.The groundwater table on the left side of the vadose zone is slowly increased to the height of h 1 = 3 m within the time t 1 = 10 days; thereafter, it is kept constant for a duration of 30 days.Afterwards, the groundwater table is further increased in height for another t 2 = 10 days, and then it is kept constant (see Figure 3).This leads to water flow in the x-direction.Figure 4 shows the change in |w LS | in the right corner of the vadose zone, which suits the change in the groundwater table on the left-hand side of the vadose zone.|w LS | increases each time with the rise of the groundwater table.When the groundwater table is constant, |w LS | remains constant, too. Figure 5 illustrates the distribution of the saturation degree s L and the pore fluid flow at different times prior to the stationary state.The water table gradually starts to increase with the increase of the groundwater table on the left side of the vadose zone.In the stationary state, the pore liquid saturation reaches its final distribution.The mesh has to be strongly refined in the zones with high gradients of effective saturation.The contaminant source was only available in the initial 25 days; it was subsequently removed.Figure 6 demonstrates the contaminant distribution within 140 days.During the first 40 days, the contaminant only moves downwards in the y-direction, without a gradient in the x-direction.The rising water table has no significant influence on the contaminant transport.Afterwards, the contaminant reaches the lower part of the studied vadose zone and begins to move in two directions, namely x and y.The maximum flux of the contaminant in the x-direction occurs in the lowest 8 m of the vadose zone, which is saturated after 78 days.The results indicate that the contaminant reached the groundwater after 76 days, in which 83% of the initial concentration entered the vadose zone.In the homogeneous horizontal contaminant transport situation, the plume must conform to an idealized wedge-like shape; however, in the vadose zone with the assumed boundary condition in this numerical example, the plume is completely deformed.When the contaminant flows over the vadose zone, there is no corresponding distortion for the wedge in the overall geometry of the system.

Conclusions
A comprehensive triphasic TPM model was developed to predict contaminant transport through the vadose zone.The developed model was applied to a numerical example, in which the contaminant source was removed after 25 days.The water level on the left boundary of the model rose from 0 to 6 m in 30 days.The results show that the contaminant has no flux in the x-direction at the highest level of 12 m.The contaminant flux increases as the contaminant moves downwards in the remaining 28 m.The concentration of the contaminant reaching the groundwater was 17% less than that of the contaminant source.The contaminant reached the ground water after 84 days.The study shows that the contamination of the vadose zone can be seen as a contaminant source for groundwater.The next step of this study could be the simulation of the filtration effect of the vadose zone and the reactive transport of contaminants through the vadose zone.The distribution of the pore gas suction is only possible if the formulation is triphasic.Moreover, the maximum gas suction can substantially affect the liquid flow situation.Despite the fact that it diminishes with time due to the gas inflow through the specimen's top surface, this effect demonstrates the most significant distinction between the triphasic and biphasic formulations.The water efflux can also demonstrate the distinction between triphasic and biphasic formulations.Furthermore, the entropy inequality is enlarged by the balance of the mass of the contaminant coupled with the Lagrange multiplier λ LC .It reads as follows: The total sum of the entropy inequality for the proposed triphasic model with a liquid contaminant, in connection with the constraints and adaptions, reads In order to evaluate the entropy inequality based on the approach used by Coleman and Noll [54], the following dependencies for the Helmholtz free energy functions are postulated: Furthermore, the capillary pressure can be defined as a difference between the pressures of the non-wetting and wetting phases using the equation below: p C = p GR − p LR , (A25) which can be reformulated by utilizing the free energy and the saturation degree of the water phase as follows: Since the change of saturation degree results in the change of the capillary pressure at the macroscopic level, p c can be defined as a function of the saturation degree, expressed as follows: The definition of pore fluid mobilities needs a relation between the liquid saturation s L and the capillary pressure p C , which can be expressed by the van Genuchten model: s L eff p C = 1 + α gen p C j gen h gen , (A28) where α gen , j gen , h gen are the material parameters, while s L eff is the effective saturation function and describes the area between the two residual saturations [55], which is defined with the equation below: Here, s L res and s G res are the residual saturation of water and air, respectively.utilizing Equation (A28), the capillary pressure can be rewritten as follows: With an assumption of s L eff ≈ s L , the Helmholtz free energy for water can be derived with The following ansatz functions are considered for the soil and contaminant, respectively:

Figure 1 .
Figure 1.Representative elementary volume (REV) of the vadose zone and its homogenization.
where C S = F T S F S and D S denote the right Cauchy-Green tensor and the the symmetric part of the deformation velocity.The momentum production term of the solid can be substituted bypS = − pL − pG − pLC .(A11)Moreover, s L L , which is required to derive ψ L L , is achieved by utilizing the definition of the saturation function with s
also be neglected.Hence, the seepage velocity of the water reads as followsn L w LS = − k L r K S µ L [grad p LR − ρ LR b].(A49)

Table 1 .
Material parameters and initial values.