Three-Dimensional Model for Bioventing: Mathematical Solution, Calibration and Validation

: Bioventing is an established technique extensively employed in the remediation of soil contaminated with petroleum hydrocarbons. In this study, the objective was to develop an improved foundational bioventing model that characterizes gas flow in vadose zones where aqueous and non-aqueous phase liquid (NAPL) are present and immobile, accounting for interphase mass transfer and first order biodegradation kinetics. By incorporating a correlation for the biodegradation rate constant, which is a function of soil properties including initial population of petroleum degrader microorganisms in soil, sand content, clay content, water content, and soil organic matter content, this model offers the ability to integrate a specific biodegradation rate constant tailored to the soil properties for each site. The governing equations were solved using the finite volume method in OpenFOAM employing the “porousMultiphaseFoam v2107” (PMF) toolbox. The equation describing gas flow in unsaturated soil was solved using a mixed pressure-saturation method, where calculated values were employed to solve the component transport equations. Calibration was done against a set of experimental data for a meso-scale reactor considering contaminant volatilization rate as the pre-calibration parameter and the mass transfer coefficient between aqueous and NAPL phase as the main calibration parameter. The calibrated model then was validated by simulating a large-scale reactor. The modelling results showed an error of 2.9% for calibrated case and 4.7% error for validation case which present the fitness to the experimental data, proving that the enhanced bioventing model holds the potential to improve predictions of bioventing and facilitate the development of efficient strategies to remediate soil contaminated with petroleum hydrocarbons.


Introduction
Bioremediation is an effective strategy for reducing concentrations of organic contaminants in environmental media such as soil, groundwater, and sediments to meet safety criteria.This process stimulates indigenous microorganisms through the addition of electron donors/terminal electron acceptors, nutrients, and growth substances, leading to the degradation of pollutants into less hazardous compounds.By-products of bioremediation include carbon dioxide, water, biomass, methane, and energy [1,2].
This eco-friendly approach capitalizes on the natural ability of microorganisms to degrade pollutants, requiring minimal effort and avoiding disruption of natural flora.In situ bioremediation eliminates health and environmental risks associated with off-site transportation of hazardous wastes.Cost-effective and sustainable, bioremediation transforms contaminants into environmentally friendly compounds, thereby eliminating future disposal responsibilities [3][4][5].Overall, it is a simple, inexpensive, and easily implemented method for cleaning up contaminated sites [6].
In bioventing (BV), as an in-situ bioremediation technique, air and nutrients are injected into polluted soil through wells to stimulate the native microorganisms.In this method, the supplied air flow rate is as low as just providing the required oxygen for the biodegradation, therefore, the dispersion and evaporation of pollutants to the atmosphere are minimized [7].The method is specifically suitable for remediation of hydrocarbons that exist in contaminated soils [8].
There are several mathematical models described in the literature that explain the movement and biotransformation process of pollutants in the subsurface.Models that only consider the behavior of gas flow are utilized for evaluating and designing BV systems.Mohr and Merz [9], employed a two-dimensional model to assess the axisymmetric air flow surrounding a vertically installed soil venting well.The primary objective was to calculate the well's radius of influence within the soil, representing the maximum distance from the well where the remediation rate remains sufficiently high to achieve a predefined cleanup goal.However, in this study there is a lack of incorporating biodegradation and the model is limited to 2D geometries.
Shikaze et al. [10] developed a two-dimensional finite element model whereby the porous matrix was represented by rectangular elements and phase partitioning was assumed to be at equilibrium.Armstrong et al. [11] developed a 2D mathematical model considering convection for the gas phase and rate limited mass transfer between phases.However, no residual NAPL phase was assumed which means the contaminant is dissolved in the aqueous phase, and the soil moisture was considered constant.The governing equations were solved applying Galerkin finite element with triangular elements and linear basis functions.Some models are limited by local equilibrium assumption (LEA) between phases which means that rate limited mass transfers of components between phases are not considered [12][13][14].McClure and Sleep [13] developed a two-dimensional model for bioventing through the application of a finite difference method.The model encompasses equilibrium interphase mass transfer within the context of three-phase (gas-organic-water) flow, considering the dispersive transport of oxygen, carbon dioxide, and organic components.The comprehensive model accounts for component consumption, oxygen utilization, microbial growth, and the resultant production of carbon dioxide during microbial activation.The Dual Monod method was employed to emulate limitations in biodegradation imposed by oxygen and substrate availability.Simulation outcomes underscored the pivotal influence of location and gas flow rates on the efficacy of the bioventing technique.While the model demonstrates thorough consideration of various bioventing-related factors, the inclusion of all three phases as mobile entities entails a considerable computational expense.Moreover, the use of the finite difference method, though effective, exhibits sensitivity to grid size, necessitating a relatively small grid for result accuracy.The finite difference method may pose limitations in addressing intricate geometries or specific boundary conditions.
Rathfelder et al. [15] presented the most sophisticated model, the Michigan Soil Vapor Extraction Remediation Model (MISER), which includes non-equilibrium interphase mass transfer and Monod kinetics for aerobic biodegradation, however MISER is limited to 1D or 2D simulations and the governing equations are solved on a structured triangular mesh incorporating finite element method which is not a mass conservative method.Additionally, the mathematical solution used in MISER is pressure-based formulation, which may cause mass balance error in the model [16].
The governing equations employed by Rathfelder et al. [15] are widely used to mathematically model SVE and bioventing.Schwarze et al. [17] conducted a modeling study on bioventing, which featured the inclusion of horizontal venting wells to elucidate the flow dynamics of pressurized air within the soil.The model deliberately omitted the consideration of biodegradation, focusing only on evaluating the physical characteristics of air movement in soil environments.The governing equations were discretized through the finite volume method.The model accounted for soil heterogeneity and variations in permeability while excluding random permeability fluctuations.
Sui et al. [18,19] proposed a two-dimensional bioventing model tailored for the remediation of trichloroethylene (TCE) and toluene contaminants in soil.The model incorporated non-equilibrium interphase mass transfer, multicomponent flow, and aerobic biodegradation, employing an upstream-weighted multiple cell balance method and operator splitting method (OS) for solution.The air injection rate emerged as a pivotal factor influencing contaminant distribution, with lower flow rates demonstrating greater efficacy in the later stages of bioventing.Notably, a substantial 95% remediation rate for total TCE was achieved.This model was limited to 1D and 2D simulations and applying finite element method.
Based on the literature, this study aims to present new features to cover some of the existing lacks and challenges.Innovations and contributions of this study are summarized as:

•
First Comprehensive bioventing mixed model using Finite Volume Method Incorporating soil properties to simulate specific sites.Employing mixed Pressure-based and Saturation-based method.
Performing simulations with any type of boundary condition, initial condition and structured or non-structured mesh.
Robustness and low computational cost.
The proposed mixed mathematical model considers the flow of the gas phase in the soil, along with the presence of immobile aqueous and NAPL phases.It also considers the transport of components in these phases and the biodegradation that occurs in the aqueous phase.Additionally, the model can be used to examine soil vapor extraction, gas flow dynamics in soil, and the effectiveness of injection and extraction wells.
One of the main challenges in implementing bioremediation is accurately predicting the contaminants degradation rate when dealing with specific site conditions.Different sites have different types of soil characterization, which makes it hard to transfer the knowledge and results from one site to another.
To address this issue and help with cost estimation for soil remediation projects, researchers in our group at the School of Engineering, University of Guelph, have been working on developing tools to estimate degradation rate-constants without the need for lengthy experiments.Khan [20] and Mosco [21] investigated BV process in meso-scale and macro-scale to simulate the field conditions in a more proper way.These studies led to a general two stage correlation between soil properties and biodegradation rate constant.
The correlated rate-constant is used in this study providing the ability to estimate the degradation rate-constant based on soil type and characterization to improve the adequacy of the closure time prediction for each specific site.This will be helpful to practitioners and consultants to make informed decisions and accurately estimate the time and cost required for remediation projects.The correlation incorporated to the proposed model is described in Section 3.4.
After developing a mathematical model, the selection of an appropriate numerical solution method is crucial to effectively address the complexities arising from equations and geometry.The advancement of bioventing simulation heavily relies on this aspect.In our previous study, Khodabakhshi & Zytner [22], different mathematical methods were compared, and it became evident that the Finite Volume Method (FVM) possesses desirable characteristics, including mass conservation, suitability for unstructured polygonal meshes, and ease of implementation for various boundary conditions.Consequently, the FVM was selected as the method of choice for this study.
To implement the simulation, OpenFOAM was specifically chosen for several of its features.OpenFOAM allows for the deployment of simulations on parallel computers and provides extensive customization options, resulting in a high degree of flexibility.Moreover, its dimensionally independent nature, coupled with automatic equation discretization, further enhances its suitability for this study.These features enable easy adjustment of parameters and equations as required, facilitating seamless integration with user-defined meshes and geometry [23].Notably, OpenFOAM supports the utilization of solvers on meshes of varying dimensions (1D, 2D, or 3D) [24], thereby accommodating a wide range of simulation scenarios in soil bioremediation.
For the numerical implementation, OpenFOAM's PMF toolbox developed by Horgue et al. [24] was improved to simulate bioventing system.This toolbox enables adaptability to various geometries and boundary conditions in three-dimensional simulations.
In a bioventing system, the primary process responsible for reducing contaminant concentration in soil is biodegradation.Due to the relatively low airflow rate, the volatilization of hydrocarbons decreases significantly compared to soil vapor extraction (SVE) [25].However, volatilization still contributes to the reduction of Non-Aqueous Phase Liquid (NAPL) concentration in soil as the second important factor.To accurately model this system, the NAPL-gas mass transfer coefficient was initially calibrated, followed by adjusting the NAPL-aqueous mass transfer coefficient using a linear function.The calibrated model was then used to simulate a large-scale reactor to validate its effectiveness.The results demonstrated that the improved mathematical model and its solution can predict the behavior of bioventing.

Conceptual Model
The soil system is modeled to include three fluid phases including a mobile gas phase, an immobile organic liquid phase (NAPL), and an immobile aqueous phase (Figure 1).The immobility of the organic liquid and aqueous phases means that convection is not considered for these phases and only biodegradation and interphase mass transfers affect their saturation in soil.
features.OpenFOAM allows for the deployment of simulations on parallel comput provides extensive customization options, resulting in a high degree of flex Moreover, its dimensionally independent nature, coupled with automatic eq discretization, further enhances its suitability for this study.These features enab adjustment of parameters and equations as required, facilitating seamless integratio user-defined meshes and geometry [23].Notably, OpenFOAM supports the utiliza solvers on meshes of varying dimensions (1D, 2D, or 3D) [24], thereby accommod wide range of simulation scenarios in soil bioremediation.
For the numerical implementation, OpenFOAM s PMF toolbox developed by H et al. [24] was improved to simulate bioventing system.This toolbox enables adap to various geometries and boundary conditions in three-dimensional simulations.
In a bioventing system, the primary process responsible for reducing conta concentration in soil is biodegradation.Due to the relatively low airflow ra volatilization of hydrocarbons decreases significantly compared to soil vapor ext (SVE) [25].However, volatilization still contributes to the reduction of Non-A Phase Liquid (NAPL) concentration in soil as the second important factor.To acc model this system, the NAPL-gas mass transfer coefficient was initially cali followed by adjusting the NAPL-aqueous mass transfer coefficient using a linear fu The calibrated model was then used to simulate a large-scale reactor to valid effectiveness.The results demonstrated that the improved mathematical model solution can predict the behavior of bioventing.

Conceptual Model
The soil system is modeled to include three fluid phases including a mob phase, an immobile organic liquid phase (NAPL), and an immobile aqueous phase 1).The immobility of the organic liquid and aqueous phases means that convectio considered for these phases and only biodegradation and interphase mass transfer their saturation in soil.The initial spatial distribution and composition of the NAPL are user-defined inputs, while the gas phase is mobile and can flow in response to applied pressure at extraction/injection wells.A component transport model is employed, where the transport and transformation of the organic liquid components are modeled.The organic liquid contaminant is considered as one component (gasoline) with the partitioning of gas and aqueous phase constituents into the organic liquid assumed negligible.The composition of the organic liquid can vary in space and time due to mass exchange into adjacent phases.
Microbes are assumed to only metabolize from the aqueous phase, with first-order kinetics incorporated into the model and the biodegradation rate constant based on a two-stage correlation developed by Khan et al. [26], which applies soil characterization into the model.
The gas phase is assumed to be comprised of air, water vapor, and the volatile components of NAPL.Mass transfer expressions are incorporated into the model to simulate rate-limited mass exchange between phases.
The prediction of transient temperature effects requires the solution of an energy balance equation, which substantially increases model complexity and computational cost.To limit model complexity, constant temperature is assumed, and all the physical and chemical properties incorporated into the model are based on 20 • C.

Phase Mass Balance Equations
The governing equations are capable of handling multi-components in three phases (gas (g), aqueous (a), NAPL (n)).Based on each simulation scenario, these components can be contaminant, Oxygen, Nitrogen, nutrient(s), water, or even a tracer.
To describe gas phase flow in soil, a general form of mass balance equation of a fluid in porous media is employed [15,18,27,28].
g represents the gas phase and c represents the components ϕ = soil porosity [-], known constant parameter S g = gas phase saturation in soil [-], unknown variable dependent on time ρ g = gas mass density [ML −3 ], known constant V g = pore velocity of the gas phase [LT −1 ], unknown dependent variable, a function of pressure ∑Eg = rate of interphase mass transfer of all components to the gas phase from other phases per unit pore volume [ML −3 T −1 ], known parameter, updated in each time step Q g = internal source/sink of the gas phase-volume of the gas phase per unit aquifer volume per unit time [T −1 ], known constant The phase velocity is described by Darcy's law as follows [15].
q g = specific discharge of the gas phase [LT −1 ] k ii = intrinsic permeability tensor of soil [L 2 ], known parameter k rg = relative permeability of the gas phase [-], known parameter P g = gas phase pressure [ML −1 T −2 ], unknown variable dependent on time and space g = gravitational acceleration vector [LT −2 ], known constant Then Equation (1) can be written as: Phase mass balance for aqueous phase (a) and NAPL phase (n) are presented based on the assumption that these phases are incompressible, immobile, and the changes in saturation of the liquid phases are only due to interphase mass transfer and biodegradation reaction.
S n = NAPL phase saturation in soil [-], unknown variable dependent on time ρ n = NAPL phase mass density [ML −3 ], known constant ∑En = rate of interphase mass transfer of all components from NAPL phase to other phases per unit pore volume [ML −3 T −1 ], known parameter, updated in each time step The impact of sorption and desorption on the bulk density of the solid phase is deemed insignificant and the bio-phase is assumed to have the same properties and volume as the aqueous phase, and as such, the solid and bio-phase do not require phase mass balances.
In Equations ( 3) and ( 5); S g , S a , S n are unknown variables which are lumped values only dependent on time.P g is a dependent variable which is a function of time and space.ϕ, ρ n , ρ a , g, k ii are constant parameters and user defined.Gas relative permeability (k rg ) is a function of the aqueous phase saturation in soil (Section 3.5).Mass transfer terms (∑E n , ∑Ea, ∑Eg) are linear functions of components concentration in each phase (Section 3.3).The net rate of biological transformation of each contaminant component (B a,c ) is also a function of concentration of contaminant component c in aqueous phase (Section 3.3).Starting with initial values of concentrations, these terms are updated in each time step after solving transport equations (Section 3.2).

Component Mass Balance Equations
To describe the fate of each component in the gas phase, a macroscopic equation that governs the transport of components in the gas phase through unsaturated soil, utilizing Darcy velocity is employed which accounts for dispersion and diffusion.Interphase mass transfer for each component is incorporated through rate-limited mass transfer kinetics [15,18,27].
C g,c = mass concentration of component c in the gas phase Component c mass balance for aqueous phase and NAPL phase are described in Equations ( 8) and ( 9) respectively.Biodegradation is assumed to occur for each contaminant component c in the aqueous phase.∂ ∂t

∂ ∂t
C n,c = mass concentration of component c in aqueous phase [ML −3 ].For each component c, an equation in each phase must be solved.

Mass Transfer
The dual resistance model is commonly used to describe non-equilibrium interphase mass transfer processes.According to this model, mass transfer is influenced by diffusion rates on each side of the interface, assuming no resistance occurs at the interface itself.Typically, the diffusional resistance on one side of the interface is the most significant, allowing the mass transfer to be characterized by an overall mass transfer coefficient [29].
The rate of organic phase volatilization is assumed to be controlled by the gas phase resistance and is described with a linear driving force expression.
C e ng,c = gas phase mass concentration of component c in equilibrium with the NAPL phase obtained from Raoult's Law [ML −3 ]  Harper et al. [30] calculated a range of volumetric mass transfer coefficient between NAPL and the gas phases with middle range soil water content as 5 h −1 to 300 h −1 .This range is converted to a mass-based range (3.14 h −1 to 189 h −1 ), using mass density of the soils and the contaminant, and was used in this work to calibrate the model prediction against accumulative volatilization experimental data.
For the dissolution of the NAPL in the aqueous phase, a linear driving force model is employed, similar to the previous scenario.In this model, it is assumed that the resistance to mass transfer is primarily present within the aqueous phase.To complete the expression, Equation (11) establishes a relationship between the equilibrium concentration of component c in the aqueous phase and the corresponding mass concentration in the NAPL phase.Harper et al. [30] suggested that instead of a constant value for mass transfer, incorporating the coefficient as a linear function of the amount of contaminant in soil helps to achieve more realistic result.This function is set as follows for NAPL-aqueous mass transfer coefficient and is used to calibrate the model against the experimental results: K an,c = mass transfer coefficient between NAPL phase and aqueous phase [ In the context of mass transfer between the aqueous and the gas phases, it is often assumed that the controlling resistance lies within the aqueous phase.This assumption holds true for organic compounds with low solubility.However, for moderately soluble compounds, the resistance in the gas phase becomes more significant [31,32].The transfer of mass between the aqueous and the gas phases can be expressed by the following equation.
K ag,c = lumped gas-aqueous mass transfer coefficient [T −1 ] C a,c = mass concentration of component c in the aqueous phase [ML −3 ] C e ag,c = aqueous phase mass concentration of component c in equilibrium with the gas phase obtained from Henry's Law [ML −3 ]  In the study conducted by Gidda et al. [33], it was proposed that the mass transfer coefficient between the aqueous and air phases for hydrocarbons falls within the range of 1 s −1 to 0.001 s −1 .The authors further noted that these values align favorably with the constant value of 0.01 s −1 suggested by Thompson and Brush [34], as well as with the reported volumetric values ranging from 0.08 s −1 to 0.009 s −1 by Harper et al. [30].Therefore, a constant value of 0.01 s −1 is incorporated in the model for aqueous-gas interphase mass transfer coefficient for the contaminant.
To address sorption in bioventing, an equation suggested by Gidda et al. [33] is employed.
Gidda et al. [33] found that a constant value was most effective for the mass transfer coefficient (k as ) from aqueous phase to soil (sorption).Therefore, k as was maintained at a constant value of 3.6 h −1 in all simulations.
Considering Gidda et al. [35] findings, it was determined that the transfer of NAPL to sorbed phase and vapor to sorbed phase did not significantly contribute to the remaining contaminant concentration.Hence, these components were excluded from the simulations.Additionally, in this study it is assumed that biophase is not a separate phase and biodegradation occurs only in the aqueous phase, therefore no mass transfer rate is considered to or from biophase.During bioventing, the mass transfer rate between all pairs of phases was regularly updated for all simulations.

Biodegradation
The model being proposed incorporates first-order kinetics since it is assumed that the degradation rate is proportional only to the concentration of the contaminant and not further limited by oxygen and/or nutrients due to abundance of oxygen and nutrients supplied during industrial bioventing.The microbial reaction is influenced by the availability of microbes to access the contaminant, which is represented in the model through the rate limited mass transfer of the contaminant into the aqueous phase.This aspect of the model is calibrated against experimental data to ensure its accuracy.
Additionally, the model employs a correlation between biodegradation rate-constant and soil properties to integrate soil characteristics into the model.A two-stage rate constant correlation derived by Khan et al. [26] is presented in Equations ( 16) and (17).
Based on the outcomes of the experiments conducted by Khan [19] and Mosco [20] for different types of soil in different reactor sizes, it is evident that a distinct alteration in degradation rates occurs within a period of 7-9 days from the start of the treatment.To facilitate comparison, Khan [19] and Mosco [20] identified Day 8 as the critical juncture where degradation rates appeared to undergo a shift, and exponential degradation curves exhibited the highest conformity.This suggests the presence of a two-stage degradation process.
In addition to the degradation patterns depicted by Khan and Mosco [20,21], a rapid increase in the microbial population was observed during the initial eight days, followed by minimal changes in the population of petroleum degrading bacteria (PDB).This observation lends further support to the concept of two-stage degradation, as also highlighted in the study by Khan et al. [22,33].In their research, the highest concentration of petroleum degrading bacteria was observed between days 7 and 8, after which the PDB population remained relatively stable for the duration of the treatment.We can consider Stage-1 as the period when favorable conditions allow indigenous microbes to become highly active due to an abundance of nutrients in the soil.Stage-2, on the other hand, represents a time when the microbial population reaches a stable level, enabling hydrocarbon degradation to occur at a realistic and balanced rate.This balanced state, which lasts for the remaining treatment days, is characterized by continuous and stable microbial functioning [26,36,37].
The correlations for the biodegradation rate constant are as follows, PDB = initial population of petroleum degrader microorganisms in soil, (CFU/g of soil) Sand = soil sand content, w% Clay = soil clay content, w% SW = soil water content, % OM = soil organic matter content, % During the initial stage, gasoline present in the soil water or weakly attached to soil particles undergoes rapid breakdown and dissolution facilitated by microorganisms.However, the second stage of degradation (from Day 8 until the experiment's conclusion) holds greater significance for the treatment as it provides valuable insights for estimating the time required for remediation.The degradation curves observed during this stage are particularly useful for developing closure times [20,21].

Gas Permeability
Gas permeability in soil can be estimated as a function of intrinsic permeability and gas relative permeability [38].
k ii = soil intrinsic permeability k rg = relative gas permeability in soil which is a function of saturation With relative gas permeability estimated using Equation (20), λ = Brooks-Corey pore size distribution index S e = effective water saturation and is calculated as follows, S a = aqueous phase saturation S r = residual water saturation.

PMF Solver Improvements
The "porousMultiphaseFoam v2107" (PMF) toolbox, initially developed by Horgue et al. [28], is a robust computational tool designed to solve the equations governing multiphase flow, employing a generalized Darcy's law, in porous media or groundwater flows using the Richards' equation.This toolbox leverages the capabilities of OpenFOAM, providing automatic discretization on three-dimensional unstructured grids and demonstrating excellent parallel efficiency.
To expand the applicability of PMF and address more complex scenarios, such as rapid water flows and solute transfers in realistic hydrological configurations with varying forcing conditions (heterogeneous infiltration and localized tracer injection), Horgue et al. [24] made significant enhancements to the toolbox.These extensions encompass coupled water flow and solute transport, enhanced numerical techniques tailored for strongly nonlinear problems, libraries, and executables for preprocessing input data, including geographical information and time-variable forcing terms, and passive or coupled scalar transport capabilities within groundwater solvers, accommodating any number of species.
With the array of solvers and tools offered by the PMF toolbox, it stands as a dependable tool for simulating fluid flow in porous media, especially when the equations involve multiphase flow and component transport.However, for simulating bioventing systems, specific modifications and advancements are necessary to encompass all aspects of bioventing.
Firstly, an essential improvement involves simulating gas flow (compressible) within porous media instead of water flow (incompressible).PMF provides a mixed formulation for Richard's equation, applicable to both saturated and unsaturated soil systems to model groundwater movement in soil.This approach is adapted and modified in this study to solve the convective equation for the gas phase in the presence of an aqueous phase and NAPL phase within soil utilized for a bioventing system.Due to low pressure difference in bioventing system (5 Pa-50 Pa), assuming that gas is an incompressible fluid seems to be reasonable [39,40].However, this model accounts for gas compression, via a gas density correction method that is used in a way, so that at the end of each time step the gas density is corrected as a function of calculated pressure distribution using an ideal gas formulation.
Secondly, the PMF toolbox currently solves the mass balance equation for only one saturation, primarily focusing on gas mass balance (Equation ( 3)).However, in a bioventing system, the presence of all three phases (gas, aqueous and NAPL) is vital.Computing the saturation of aqueous phase is important because it is assumed that biodegradation occurs only in aqueous phase, and tracking the changes in saturation of NAPL phase leads to finding the closure time.To incorporate these equations, solutions for two Ordinary Differential Equations (ODEs) addressing Equation (4) for the aqueous phase and Equation (5) for the NAPL phase are added to the main code of the solvers.
Thirdly, a critical enhancement involves integrating interphase mass transfer terms into the main solver.This integration is achieved by incorporating source/sink terms into each equation.These terms are updated in each time step, reflecting the concentration values, and are utilized as constant values in the subsequent time step, representing mass transfer terms.

Nonlinear Gas Mass Balance Equation Solution
In order to address the equations governing fluid flow and component transport, the initial gas pressure distribution and Darcy velocity must first be determined.According to existing literature, given that the equation aims to be solved for single-phase flow (gas) in the presence of an aqueous phase, several mathematical approaches for the mass balance equation (Equation (3)) exist.In a mass balance equation for each mobile phase, the phase pressure and phase saturation are unknown variables and by choosing one or both as the primary variable, different approaches are described.
In the pressure-based method, fluid pressures are selected as the primary dependent variable and the flow equations are solved simultaneously for the aqueous phase pressure (P a ) and the gas phase pressure (P g ).Saturations are subsequently updated from capillary pressure-saturation relations.The approach offers broad applicability for variably saturated flow and heterogeneous soils.This is primarily due to the monotonous nature of hydraulic head across the variably saturated domain, as well as its continuity across the interface of heterogeneous soils.Rathfelder et al. [15] implemented this method to solve the flow equations for two-phase flow (gas and aqueous) in soil to simulate SVE and BV.Nevertheless, the non-linear capillary pressure term may introduce mass balance errors in the numerical model, as highlighted by Celia et al. [16].
The saturation-based method, where phases saturations are selected as primary variables, proves highly efficient in simulating unsaturated flow in homogeneous soil, as the resulting numerical model inherently maintains mass conservation and exhibits lower non-linearity compared to the pressure-based model under dry soil conditions.However, this method lacks the capability to handle locally saturated areas, a feature present in the pressure-based method, as noted by Zha et al. [41] and Horgue et al. [24].
In light of these considerations, Celia et al. [16] proposed a mixed-form solution that combines the advantages of both pressure-based and saturation-based models.The mixed formulation ensures mass conservation akin to the saturation-based form while simultaneously accommodating locally saturated areas as observed in the pressure-based form.In this study, the mixed formulation is employed to solve the mass balance equation for the gas phase.
Furthermore, Horgue et al. [28] created a general open-source toolbox using Open-FOAM for generic two-phase flow, incorporating widely employed soil water retention models for relative permeability and capillary pressure.This toolbox, known as PMF, was later expanded to accommodate Richards' equation by Horgue et al. [42], utilizing pre-existing libraries and functions, thereby delivering an effective solution specifically tailored for air-water flows.Additionally, new solvers were developed and added to PMF by Horgue et al. [24].In the present work these solvers from PMF were improved to adopt for bioventing system.

•
Stationary unsaturated flow in porous media.

•
Passive multi-scalar transport solvers on pre-computed flow fields in 3D.

•
Coupling 3D flow modeling with multi-scalar tracer transport.
In this section the solution method for the nonlinear convection equation is explained.Equation ( 3) can be written as follow, where, In situations where interphase mass transfer to and from the gas phase is deemed low and exerts negligible influence on the gas flow dynamics, the mass transfer component on the right-hand side of the equation may be disregarded.Thus, the density of gas is defined with an ideal gas model.
Considering the definition of phases saturations in soil, and neglecting the effect of changes in NAPL saturation on gas flow due to low NAPL saturation (~0.01), Now, using capillary pressure definition saturation and pressure can be related as, C p is calculated using capillary models by Van Genuchten [43] which is added to PMF toolbox by Horgue et al. [28].
To reduce the gas mass balance equation to a sequence of linear equations, it is required to express the implicit gas saturation (S g ) at the iteration (i + 1) as a function of gas pressure (P g ) by utilizing capillary pressure.
i + 1 represents the next iteration while i shows the parameters in the current iteration.
S g represents the gas saturation at the previous time step.This allows the mass conservation equation to be expressed as a function of the solved variable P g i+1 along with the current values of P g i and S g i .
Density of the gas, ρ g is assumed to be constant in each time step.k rg is also a function of aqueous phase saturation (S a ), however, since it is crucial to maintain water content in an optimum range during bioventing process [20,21], it is reasonable to assume that relative permeability (k rg ) is constant.Now the equation is reduced to: Equation ( 31) is solved iteratively using Picard's algorithm, where the matrix A i contains accumulative terms, P g represents the vector of gas pressure at the previous time step, M i consists of terms related to pressure gradients, L i involves terms related to gravity, and G i represents the vector of source term and mass transfer rates which are a function of concentrations calculated in previous time step (in the first time step initial values of concentrations are used to calculate the mass transfer rates).During each iteration, the matrix coefficients are updated using the pressure vector P g i at that iteration i.The residual of the equation at iteration i is denoted as R i .
The iteration is terminated and assumed converged if, and r tol,Pic is a user defined value.Picard's iteration method is robust and has a wide range of convergence for low flow rates.However, to enhance the convergence rate and reduce computational costs for highly nonlinear problems, Horgue et al. [24] suggested incorporating Newton's method sequentially into the model.Due to the complex behavior of the Van Genuchten function relating the gas saturation and pressure, especially near the saturated point where the capillary function derivative tends to infinity as saturation approaches 1, the process of Picard's algorithm encounters difficulty in achieving convergence.This situation necessitates numerous iterations and small-time steps, which leads to a considerable computational burden.Consequently, simulating this configuration within a reasonable timeframe becomes impractical, with execution times varying from under 5 min when employing the newly integrated Newton's algorithm to several hours using Picard's algorithm on a standard CPU.
This solution dilemma is overcome by initially solving the equation using Picard's method until the residual |R i | reaches a tolerance value, r tol,Pic .Subsequently, the solution is refined using Newton's method until the residual |R i | is below a desired tolerance value, r tol,New .This approach leverages the broad convergence range of Picard's algorithm to approach the solution roughly, followed by the quadratic convergence of Newton's method to converge more accurately.
The determination of suitable convergence criteria for the mass balance equation and the component transport equation has been a subject of thorough investigation within the modeling literature.Numerous studies have delved into the intricate details of each model parameters to establish precise convergence residual values.Notably, Rathfelder et al. [15] undertook an in-depth exploration of numerical errors and recommended a convergence tolerance of 10 −4 to 10 −3 for the pressure equation, and for the component transport equation a magnitude of 10 −10 .Williams and Miller [44] investigated the convergence of different approaches for solution of Richard's equation, proposing an optimal range of 10 −3 to 10 −2 for the pressure equation, adapting it to varying soil properties and model parameters.
Sin et al. [45] contributed significantly by investigating the intricate implementation of a compressible multicomponent two-phase flow utilizing a reactive transport simulator.Through their research, they advised a convergence range of 10 −8 to 10 −6 for the gas flow equation, 10 −5 for the reactive transport equation, and a more stringent 10 −12 for complex geochemical reactions, underscoring the multifaceted nature of their simulations.
Drawing upon this foundation, Horgue et al. [24] incorporated a convergence residual of 10 −10 for the validation cases within the PMF toolbox.In the context of our present study, we adopt the convergence residual value of 10 −10 for our simulations concerning the gas flow in soil and the component transport equation, aligning with the convergence standards set forth by these influential prior works.

Component Transport Equation Solution
The solution of the component transport equations is considered in two different scenarios in the model.In the first method, mass balance equations and transport equations are coupled and solved at the same time step.This method is capable of investigating changes in the aqueous phase and the gas phase saturations, considering transient gas flow and incorporating the effect of components interphase mass transfer on gas flow in a simulation.However, this method comes with relatively high computational cost.To do this type of simulation the following solver of the PMF toolbox is improved to be incorporated for bioventing model: • "groundwaterTransportFoam" solver; coupling 3D unsaturated flow modeling with multi-scalar component transports.
In some cases, one needs to do a simple simulation to investigate component transport and mass transfer between phases and consider a steady gas flow which does not change due to presence of other phases.Therefore, another option is added to this model to reduce the computational cost with acceptable adequacy.Users can run the simulation for the gas phase steady flow in porous media using one of the following solvers depending on the simulation conditions: • "darcyFoam" solver; Solving Darcy's equation for a single incompressible phase flow.
• "steadyGroundwaterFoam" solver; Stationary unsaturated flow in porous media in presence of second phase without considering interphase mass transfer.
The pressure distribution and velocity distribution results were then used to solve the transport equations [22].

Equation Discretization and Time Step Management
The model has been created within the OpenFOAM environment, which serves as a development framework for solving equations utilizing the finite volume method (FVM).Within OpenFOAM, there are built-in functions that automatically discretize common operators (∂/∂t, ∇, ∆) on unstructured three-dimensional meshes using a wide range of numerical schemes.Which means in each simulation it is possible to adjust the schemes to solve the equations based on the conditions of a specific scenario [23,24].
In this study, a conventional second-order centered scheme is employed for spatial discretization of convective and diffusive terms.
The calculation of the time step is performed to ensure that the relative time truncation error, estimated using the Taylor series, remains below a predefined error value (ε error ) set by the user.For the mass balance equation, backward Euler scheme (Equation ( 34)) is adapted to make the implementation of non-linear Picard and Newton's methods easier.
A second order back-ward scheme (Equation ( 35)) is used for time discretization of component transport equation, where χ max is the value of the variable at the cell giving the maximum derivative.Maximum and minimum time steps are user defined providing the option to manage the computational cost based on model parameters and mesh size.In this study, the minimum time step value is set at 20 s and the maximum value is 3600 s.The total simulated period is 30-40 days and higher time steps does not affect the computational cost significantly.

Bioventin Simulation Methodology
Considering the proposed mathematical model and its capabilities, in this section the model is utilized to simulate a specific scenario of bioventing based on the experimental works of our research group to present the calibration and validation procedure of the model against the experimental data.There are some capabilities of the mathematical model such as soil heterogeneity, incorporating Monod type kinetics, and multicomponent transport, which are not applicable in this scenario and will be investigated in future works.
In the research conducted by Khan et al. [26,36], bioventing simulations were performed using a meso-scale reactor containing 4 kg of soil based on its dry weight.Two distinct soil types were gathered from two locations, namely the Elora research station in the County of Wellington, Ontario, Canada (Elora soil) and the Delhi research farm in Delhi, Ontario, Canada (Delhi Soil).The difference between the two soils stems from their respective compositions: Delhi soil is characterized by a loamy sand composition, while the Elora soil is composed of silt loam.The soil was homogeneously contaminated with synthetic gasoline, and five distinct soil samples were examined.To enhance the effectiveness of the process, water and nutrients were added to each soil sample.Initial tests utilizing the meso-scale system demonstrated the suitability of the experimental setup for conducting reliable bioventing treatments.The parameters monitored during the treatments were the soil's total petroleum hydrocarbon (TPH) concentration, air flow rates, loss of gasoline through volatilization, and the concentration of oxygen in the off gases.
To create a model of this reactor, due to the experiment's conditions and data, the transportation of gasoline (pure NAPL phase) as a single component during the bioventing process is considered.This simulation scenario addresses the movement of the contaminant from the non-aqueous phase liquid (NAPL) to both the aqueous and the gas phases.Due to the method used in the experiments to prepare the contaminated soil, the soil system is considered homogeneous (contaminant concentration in soil, soil water content and gas relative permeability).
To ensure optimal bioactivity of microorganisms, it is crucial to maintain an adequate range for water content in the soil throughout the bioventing process.To address this, we introduce a source term into the mass balance equation of the aqueous phase.The specific value for this term is determined as reported by Mosco [21] as 253 mL water/m 3 soil per day.Additionally, it is assumed that nutrients are abundantly available, as the soil used in the experiments has been amended with nutrients.
In this scenario, the meso-scale reactor is a cylinder with dimensions of 15.4 cm × 17.8 cm (inner diameter × height, respectively).The inlet of the reactor is positioned along the outer wall, while gas flow occurs radially through the cylindrical reactor.Additionally, a venting well is situated at the center of the reactor with a diameter of 2 cm (r w = 1 cm).One-third of the well, starting from the top of the reactor, is sealed, and considered a "no flux boundary condition".The remaining portion of the well serves as the outlet.The two walls located at the top and bottom of the reactor are assumed to have "no flux boundary conditions" as well.Figure 2 shows the geometry and the boundary conditions for this simulation scenario.
In this study, ANSYS Fluent software was utilized to create the geometric structure and mesh for our investigation.A triangular mesh was generated using the auto-meshing tool feature.Strategically a finer mesh was established along the boundaries and a coarser mesh within the main computational domain.
Building upon the insights of Roman et al. [46], who proposed a guideline for simulating the movement of two-phase fluids in porous media, the goal was to represent each pore with a minimum of ten grid cells.To apply this concept, initially the minimum mesh size was set to 2.5 µm, assuming an average soil pore size of 25 µm.
In order to validate the accuracy of our simulations and ensure that the mesh size did not compromise the results, a mesh convergence analysis was performed.By varying the mesh sizes, the predicted pressure distribution was evaluated within the system.The findings indicated that a grid size of 0.25 mm yielded consistent results for gas pressure distribution in the reactor.Consequently, this mesh size was adopted, which consisted of approximately 7.23 × 10 5 cells.To facilitate compatibility between different software platforms, the mesh geometry from ANSYS Fluent's Mesh.stl format (in ASCII) was converted to OpenFOAM's Mesh.msh format using the "fluent3DMeshToFoam" command.In this study, ANSYS Fluent software was utilized to create the geometric structure and mesh for our investigation.A triangular mesh was generated using the auto-meshing tool feature.Strategically a finer mesh was established along the boundaries and a coarser mesh within the main computational domain.
Building upon the insights of Roman et al. [46], who proposed a guideline for simulating the movement of two-phase fluids in porous media, the goal was to represent each pore with a minimum of ten grid cells.To apply this concept, initially the minimum mesh size was set to 2.5 µm, assuming an average soil pore size of 25 µm.
In order to validate the accuracy of our simulations and ensure that the mesh size did not compromise the results, a mesh convergence analysis was performed.By varying the mesh sizes, the predicted pressure distribution was evaluated within the system.The findings indicated that a grid size of 0.25 mm yielded consistent results for gas pressure distribution in the reactor.Consequently, this mesh size was adopted, which consisted of approximately 7.23 × 10 5 cells.To facilitate compatibility between different software platforms, the mesh geometry from ANSYS Fluent s Mesh.stl format (in ASCII) was converted to OpenFOAM s Mesh.msh format using the "fluent3DMeshToFoam" command.
Elora and Delhi soil properties used in these simulations are extracted from [20,21] and presented in Table 1.Soil properties are used to calculate the biodegradation rate constant.
Table 1.Delhi and Elora soil properties used by Khan [20] and Mosco [21].Elora and Delhi soil properties used in these simulations are extracted from [20,21] and presented in Table 1.Soil properties are used to calculate the biodegradation rate constant.Constant parameters and initial values for the calibrating simulations are presented in Table 2.

Model Calibration
Determining the adjustable parameters for different modeling scenarios to improve the best fit does not follow strict guidelines [48].However, a general calibration principle involves modifying parameters that exhibit significant uncertainty and exert substantial influence on the sensitivity of the numerical solutions obtained from the constructed models.This adjustment process encompasses parameters that cannot be directly measured, but are acquired through model calibration, such as mass transfer coefficients.On the other hand, non-adjustable parameters integrated into the mathematical model are considered fixed values because they are measurable, have relatively low or moderate uncertainty, or have minimal impact on the results.
In terms of adjustable parameters for SVE and BV models, certain factors like mass transfer coefficients, intrinsic and relative permeability, dispersivity, porosity of soil, bulk density, and moisture content can be selected.Researchers Vogele et al. [49] and Barnes and White [50] adopted an approach where they calibrated their SVE models by adjusting the initial concentration and intrinsic permeability, focusing on these two parameters due to their high level of uncertainty as determined through evaluation.Harper et al. [30] used a linear function of NAPL-gas mass transfer equation to calibrate their model.Zhao [47] presented a complex SVE mathematical model and used exponential functions to describe mass transfer coefficients and adjusted them to calibrate the model.
In a bioventing process, biodegradation reaction is the dominant phenomena, and it is reasonable to assume that parameters affecting the reaction control the effectiveness and closure time of bioventing.In this model a correlation is incorporated to obtain the biodegradation rate constant as a function of soil properties, and it differs from soil to soil, so it is not acceptable to choose the calculated rate constant as a calibrating parameter.The other parameters which can be potentially adjusted to calibrate the model are, mass transfer coefficients.The assumption of biodegradation occurring only in the aqueous phase, makes it clear that the mass transfer between NAPL and aqueous phases is the main controlling parameter to predict the reduction of contaminant in soil [20,51].It is also important to note that NAPL volatilization is the other parameter which affects the closure time.However, both Khan and Mosco [20,21] showed that due to low gas flowrate in bioventing, volatilization of the contaminant reduces dramatically.
Importing a constant average value of reported mass transfer coefficients, the results showed a reasonable trend, but the slope of the model prediction was far from experimental data which means the closure time prediction is not acceptable.Therefore, it was crucial to calibrate the model in a physically reasonable way.
In the first step, a pre-calibration was performed by choosing NAPL-gas mass transfer coefficient as the pre-calibrating parameter to adequately predict the amount of volatilized hydrocarbon.Ten simulations were performed using 10 different values of NAPL-gas mass transfer coefficient and the results were compared to the experimental data from [20].The value of 4.3 h −1 was then chosen as the mass transfer coefficient between NAPL and the gas phases.Figure 6 shows the pre-calibrated model results compared to the experimental data for volatilized gasoline.After setting the volatilization rate, the maximum and minimum values for mass transfer coefficient of gasoline from NAPL phase to aqueous phase were set as the boundaries of calibration simulations (3.6 h −1 -36 h −1 ) [35].Ten simulations were performed, and the results were compared against the experimental data.The results still did not show a good fit in slope.After setting the volatilization rate, the maximum and minimum values for mass transfer coefficient of gasoline from NAPL phase to aqueous phase were set as the boundaries of calibration simulations (3.6 h −1 -36 h −1 ) [35].Ten simulations were performed, and the results were compared against the experimental data.The results still did not show a good fit in slope.
To address the observed error, it is reasonable to define the NAPL-aqueous mass transfer coefficient as a function of the NAPL and water content in the soil.However, in the case of bioventing, maintaining a constant water content is crucial for optimal treatment.Therefore, it is a viable choice to calibrate the model by defining the mass transfer coefficient as a function of the NAPL amount in the soil.Gidda et al. [35] used an exponential function to estimate the aqueous-air and NAPL-air mass transfer constants, as they observed a tailing effect in soil vapor extraction (SVE), leading to a significant decrease in NAPL volatilization.However, in bioventing, there is limited tailing effect, and the use of exponential terms seems unnecessary.Hence, in this study, another approach suggested by Harper et al. [30] is adopted.Harper et al. [30] proposed incorporating the mass transfer coefficient as a linear function of the contaminant amount instead of a constant value, which helps achieve more realistic results.This function is represented by Equation (12).
As previously stated, the second stage of biodegradation corresponds to a period in which the microbial population stabilizes, facilitating hydrocarbon degradation at a sustainable and steady pace.This state of equilibrium persists throughout the remaining days of the treatment, with microbial activity exhibiting steady functioning [20,21].Therefore, the data for the second stage were used to calibrate the model.
It is expected to have the maximum value of k na,c at the start of the treatment and a decrease in the value as the amount of contaminant decreases.First, the value for adjusting parameter b is set to 3.6 h −1 and a range of 25 h −1 to 50 h −1 was set for adjusting parameter a and various simulations were performed (Table 3).A log transformation was performed on the experimental data and the results from the simulation, and normalized sum of the squared relative deviations (NSSRD) values were calculated using following formula.
where N is the number of sampled points, i is i-th sampled point for experimental and predicted data.C exp,i is the concentration of contaminant in soil from experimental data and C mod,i is the concentration of contaminant in soil from model prediction results.
In the second step, the value of a was set to 45 h −1 and the same calibration process was performed to adjust the best fit for b, and the value of 2.8 h −1 was then chosen for b.Finally, Equation (38) was incorporated into the model to calculate NAPL-aqueous mass transfer coefficient in each time step.
The calibrated model was then employed to simulate the meso-scale reactor.Since the experimental data provided the batch concentration of the contaminant in soil, the saturation of contaminant S n was converted to mg of contaminant per kg of soil and summed with the accumulated amount of the sorbed contaminant on soil to capture the total amount of the contaminant in soil in each time step.Results are presented in Figure 7.The Root Mean Squared Error (RMSE) serves as a metric for quantifying the goodness of fit between model outputs and observed values.
The RMSE calculated for calibrated model is 2.9% which shows the accuracy of model prediction.

Model Validation
To ensure the model s effectiveness, it is important to test its performance through various bioventing scenarios with different conditions.This helps verify that the model s predictions remain accurate, regardless of the experimental setup.A different set of data obtained by Mosco [21] from a large-scale 80 kg reactor were used.The geometry of the rector, the boundary conditions and initial values are the same as previous simulations.Size of the reactor is changed, and air flow rate is set to 8.9 mL/min as reported by Mosco [21].
The results presented in Figure 8.The model prediction of the treatment process shows a good fit to the experimental data.The Root Mean Squared Error (RMSE) serves as a metric for quantifying the goodness of fit between model outputs and observed values.
The RMSE calculated for calibrated model is 2.9% which shows the accuracy of model prediction.

Model Validation
To ensure the model's effectiveness, it is important to test its performance through various bioventing scenarios with different conditions.This helps verify that the model's predictions remain accurate, regardless of the experimental setup.A different set of data obtained by Mosco [21] from a large-scale 80 kg reactor were used.The geometry of the rector, the boundary conditions and initial values are the same as previous simulations.Size of the reactor is changed, and air flow rate is set to 8.9 mL/min as reported by Mosco [21].
The results presented in Figure 8.The model prediction of the treatment process shows a good fit to the experimental data.Another RMSE calculation is done to present the fitness of validation case and the result shows that the error is 4.7% which is an acceptable error for the model prediction against experimental data.
Conducting simulations utilizing different geometries and environmental conditions showed the robustness and predictive accuracy of the model and associated software package in determining closure times efficiently and accurately.This model can effectively investigate and simulate bioventing and soil vapor extraction (SVE) scenarios, encompassing remediation of sites with distinct boundary conditions, and remediation involving both extraction and injection processes.Furthermore, the model considers intricate factors such as soil characteristics, making it a valuable tool for comprehensive soil remediation analyses and planning.

Conclusions
Through a thorough investigation of different models and their underlying assumptions, an enhanced mixed mathematical model was specifically designed for industrial bioventing in vadose zones.
One of the significant innovations of the model lies in its comprehensive incorporation of all relevant aspects of bioventing.By considering fluid flow in porous media in the presence of three phases, biodegradation kinetics, component transport, and soil characteristics, the model provides a holistic representation of the bioventing process.The model is advanced as it incorporates soil properties including soil organic matter content, sand content, clay content, soil water content, and initial population of petroleum Another RMSE calculation is done to present the fitness of validation case and the result shows that the error is 4.7% which is an acceptable error for the model prediction against experimental data.
Conducting simulations utilizing different geometries and environmental conditions showed the robustness and predictive accuracy of the model and associated software package in determining closure times efficiently and accurately.This model can effectively investigate and simulate bioventing and soil vapor extraction (SVE) scenarios, encompassing remediation of sites with distinct boundary conditions, and remediation involving both extraction and injection processes.Furthermore, the model considers intricate factors such as soil characteristics, making it a valuable tool for comprehensive soil remediation analyses and planning.

Conclusions
Through a thorough investigation of different models and their underlying assumptions, an enhanced mixed mathematical model was specifically designed for industrial bioventing in vadose zones.
One of the significant innovations of the model lies in its comprehensive incorporation of all relevant aspects of bioventing.By considering fluid flow in porous media in the presence of three phases, biodegradation kinetics, component transport, and soil characteristics, the model provides a holistic representation of the bioventing process.The model is advanced as it incorporates soil properties including soil organic matter content, sand content, clay content, soil water content, and initial population of petroleum degrader microorganisms in soil, through a biodegradation rate constant correlation that accurately reflects the influence of different soil types.This feature adds a layer of precision to the simulation process, crucial for tailoring remediation strategies to specific site conditions.The model serves as a reliable tool for forecasting the closure time of unsaturated bioventing, a question many clients, consultants and regulators have.The model's versatility is another noteworthy aspect.The mathematical model built in OpenFOAM is dimension independent and enables simulations in 1D, 2D, and 3D based on the imported mesh by the user, catering to the diverse needs of researchers and practitioners in the soil remediation domain.
The model excels in robustness and computational efficiency, making it accessible for simulating simple geometries on personal computers, thus minimizing the computational burden typically associated with complex simulations.Moreover, it supports a diverse array of geometries and boundary conditions, empowering researchers to execute complex simulations on parallel computing cores and cloud-based systems.
Lastly, the model's ability to simulate both Soil Vapor Extraction (SVE) and bioventing is a significant advancement, providing researchers and practitioners a comprehensive tool to evaluate and compare these soil remediation techniques, ultimately contributing to more informed decision-making in the remediation process.
To solve the nonlinear convection equation, a mixed-form solution is utilized, harnessing the advantages of both pressure-based and saturation-based models.The calculated saturation and pressure are then employed to solve the component transport equations.
For model calibration, experimental data from a meso-scale (4 kg) reactor is used.Initially, a pre-calibration is performed by considering the mass transfer coefficient between the NAPL phase and the gas phase as the pre-calibrating parameter.The predicted volatilization results are compared against volatilization experimental data to determine the best fit for the NAPL-gas mass transfer coefficient.
Once the best fit is determined, a linear equation is employed to describe the changes in the NAPL-aqueous mass transfer coefficient during the treatment process, which is a function of the amount of contaminant in the soil.The calibrating parameter's minimum and maximum values are selected based on literature, and various simulations are conducted to find the optimal fit and calibrate the model.
To validate the effectiveness of the calibrated model, another simulation is performed using initial conditions and values for a large-scale (80 kg) reactor to predict the closure time.The prediction results are then compared to the experimental data and a satisfactory fit is achieved, demonstrating the accuracy and efficacy of the proposed model and toolbox in adequately predicting bioventing closure time in the vadose zone.
In future studies, the proposed mathematical model will be employed to explore additional capabilities, such as incorporating soil heterogeneity by considering variable gas relative permeability and/or soil intrinsic permeability and layered soil, incorporating different kinetics such as Monod type kinetics, and investigating simulation scenarios considering multicomponent transport.

Figure 1 .
Figure 1.Conceptual model of bioventing system.The initial spatial distribution and composition of the NAPL are user-defined while the gas phase is mobile and can flow in response to applied press extraction/injection wells.A component transport model is employed, where the tra
S a = aqueous phase saturation in soil [-], unknown variable dependent on time ρ a = aqueous phase mass density [ML −3 ], known constant ∑Ea = rate of interphase mass transfer of all components to aqueous phase from other phases per unit pore volume [ML −3 T −1 ].Known parameter, updated in each time step B a,c = net rate of biological transformation of component c in aqueous phase per unit aquifer volume [ML −3 T −1 ], known parameter, updated in each time step ρ n ∂ ∂t C a,c = mass concentration of component c in aqueous phase [ML −3 ] B a,c = net mass rate of biodegradation of component c in aqueous phase per unit aquifer volume [ML −3 T −1 ]

K
na,c = lumped aqueous-NAPL mass transfer coefficient [T −1 ] C a,c = mass concentration of component c in the aqueous phase [ML −3 ] C e na,c = aqueous phase mass concentration of component c in equilibrium with the NAPL phase, equals to solubility limit of the component c in water [ML −3 ]
k ii = intrinsic permeability of soil, [L 2 ] Q g = flow rate of venting well, [L 3 T −1 ] P w = venting absolute pressures on the venting well, [ML −1 T −2 ] P atm = atmospheric pressure, [ML −1 T −2 ] r w = radius of venting well, [L] r I = radius of influence of venting well, [L] H = height of the extraction well through contamination region, [L] Using Equation (36), P w -P atm was calculated to be 35 Pa for this particular scenario.Figures 3-5 present the pressure distribution, streamlines and velocity distribution through the reactor, respectively.Math.Comput.Appl.2024, 29, x FOR PEER REVIEW 18 of 26

Figure 3 .
Figure 3. Pressure distribution in the reactor in 3D.Figure 3. Pressure distribution in the reactor in 3D.

Figure 3 .
Figure 3. Pressure distribution in the reactor in 3D.Figure 3. Pressure distribution in the reactor in 3D.

Figure 3 .
Figure 3. Pressure distribution in the reactor in 3D.

Figure 5 .
Figure 5.The gas phase velocity distribution in the reactor in 3D.

Figure 3 .
Figure 3. Pressure distribution in the reactor in 3D.

Figure 5 .
Figure 5.The gas phase velocity distribution in the reactor in 3D.Figure 5.The gas phase velocity distribution in the reactor in 3D.

Figure 5 .
Figure 5.The gas phase velocity distribution in the reactor in 3D.Figure 5.The gas phase velocity distribution in the reactor in 3D.

26 Figure 6 .
Figure 6.Volatilization mass rate of contaminant in different time intervals.

Figure 6 .
Figure 6.Volatilization mass rate of contaminant in different time intervals.

26 Figure 7 .
Figure 7.Comparison of calibrated model results against experimental data.

Figure 7 .
Figure 7.Comparison of calibrated model results against experimental data.

Figure 8 .
Figure 8. Concentration of contaminant in soil from experimental data and model prediction.

Figure 8 .
Figure 8. Concentration of contaminant in soil from experimental data and model prediction.
* PDB is initial population of petroleum-degrading bacteria, (log CFU/g soil), where CFU is colony-forming units.**WHC is water holding capacity of soil.

Table 3 .
Variations of NSSRD Values during Calibration.