Development of a Model for Predicting the Volatilization Flux from Unsaturated Soil Contaminated by Volatile Chemical Substances

This work developed a model for predicting the volatilization flux from the unsaturated soil contaminated by volatile chemical substances (VCSs) such as mercury and benzene. The model considers a series of phenomena under the unsaturated condition such as multi-phase flow consisting of a non-aqueous phase liquid, water, and gases together with the permeation of rainfall into the surface soil, the volatilization/condensation of VCSs, and the adsorption/desorption of VCSs. On this basis, this work clarified a mechanism for the generation of a volatilization flux at the ground surface. In addition, the effects of various transport phenomena in the surface soil on the magnitude and seasonal changes in this flux due to variations in weather factors such as rainfall level, temperature, and air pressure were quantitatively evaluated. This newly developed prediction model can be utilized to estimate dynamic variations in the flux under real-environmental conditions.


C s
The amount of adsorption of benzene to the soil particles (mg/kg) C sat,k Solubility into water phase of component k (mg/L) C w Dissolved concentration of benzene in water phase (mg/L) D A Average grain diameter (m) D A0 Average grain diameter of Toyoura sand as a standard value (m) d pore Principal pore diameter (m) Relative permeability to gas phase (dimensionless) k rg,gn Relative permeability to NAPL phase in gas-liquid two phase system (dimensionless) k 0 rg Relative permeability to gas phase, corresponding to residual liquid saturation S lr (dimensionless) k rl Relative permeability to liquid phase (dimensionless) k 0 rl Relative permeability to liquid phase, corresponding to residual gas saturation S gr (dimensionless) k rn Relative permeability to NAPL phase (dimensionless) k rn,gn Relative permeability to NAPL phase in gas-NAPL two-phase system (dimensionless) k rn,nw Relative permeability to NAPL phase in NAPLwater two-phase system (dimensionless) k 0 rn Relative permeability to NAPL phase, corresponding to irreducible water saturation S wi (dimensionless) k rw Relative permeability to water (dimensionless) k rw,nw Relative permeability to water phase in NAPLwater two-phase system (dimensionless) k 0 rw Relative permeability to water phase, corresponding to residual NAPL saturation S nr (dimensionless) k v Rate constant necessary to reach gas-liquid equilibrium (1/Pa/s) N gw Reduction order of P c,gw (dimensionless) N mol,k Total number of moles of component k (kmol) N nw Reduction order of P c,nw (dimensionless) P System pressure (Pa) P c,gn Capillary pressure operating between gas and NAPL phases (Pa) P c,gw Capillary pressure operating between gas and water phases (Pa) P c,nw Capillary pressure operating between NAPL and water phases (Pa) P d,gw Threshold pressure in gas-liquid two-phase system (Pa) P d,nw Threshold pressure in NAPL-water two-phase system (Pa) P g Pressure of gas phase (Pa) P n Pressure of NAPL phase (Pa) P sat,k Saturated vapor pressure of component k ( Irreducible water saturation obtained for Toyoura sand as a standard value (dimensionless) T Temperature (K) Surface tension in gas-water system (N/m) Φ g Flow potential of gas phase (Pa) Φ n Flow potential of NAPL phase (Pa) Φ w Flow potential of water phase (Pa) Porosity (dimensionless)

Introduction
In recent years, there have been numerous reports concerning the contamination of soil and groundwater derived from volatile chemical substances (VCSs) such as mercury (Hg) and benzene (C 6 H 6 ) as a result of leakage from factories and illegal dumping. Figure 1 presents a diagram summarizing the contamination of soil and groundwater by VCSs. These compounds tend to permeate into the pore space in the surface soil (in the region defined as the unsaturated zone) as a result of advection and dispersion and can subsequently transfer to either groundwater or the ambient atmosphere with the changes of mass and composition based on several phenomena. These include the elution of undiluted VCSs into pore water and precipitation from the aqueous phase, volatilization into gas phase and condensation into liquid phase, and partitioning between pore water and soil particle as the solid phase due to adsorption/desorption. Once these VCSs are discharged into the environment, a health risk on human body due to inhalation or intake might occur. Therefore, it is very important to understand the transport phenomena of these contaminants in soil and groundwater in order to minimize the spread of contamination and apply the appropriate countermeasures for remediation [1]. In particular, the accurate prediction of volatilization flux from the ground surface based on detailed mathematical modeling of the advection and dispersion behaviors of VCSs under unsaturated condition is necessary for the quantitative evaluation of the risks derived from these compounds.
There have been many numerical analyses assessing the contamination of soil and groundwater by VCSs. Such studies have examined the diffusion through soil of volatilized VCSs under unsaturated conditions. These works can be classified into those that tried the optimization of a series of parameters for the evaluation of the advection and dispersion of volatilized VCSs via comparisons with experimental data [2][3][4][5][6], and those that attempted to predict the distribution of volatilized VCSs targeting the soil gas extraction as remediation process and risk assessment [7][8][9][10][11][12][13][14][15][16][17][18]. Several studies also investigated the behavior of volatized VCSs near the ground surface. Briggs and Gustin [11] found that the evaporation stage/rate of mercury depended on soil characteristics such as soil particles and water content, and developed a conceptual prediction model for mercury volatilization flux. Their work indicated that volatilization flux was affected by the mass flow of water phase containing VCSs to the ground surface. Abreu and Johnson [12] developed a soil vapor intrusion model for the quantification of the attenuation coefficient (that is, indoor air concentration/source vapor concentration) based on the gas phase flow under unsaturated conditions and partitioning between gas-liquid/liquid-solid phases. The resulting mathematical model assumed that the volatilized VCSs migrate indoors through foundation cracks in buildings. Abreu and Johnson also examined the effects of vapor source-building separation, building construction and biodegradation on the change of indoor air concentration of benzene. Using Abreu and Johnson's model, Ma et al. [16] considered the Fig. 1 A schematic illustration of soil and groundwater contamination by volatile chemical substances (VCSs) effect of source concentration and soil type on the indoor air concentration of methyl tert-butyl ether (MTBE), tert-amyl methyl ether (TAME), and ethyl tert-butyl ether (ETBE) in gasoline leaking from an underground storage tank, using the odor thresholds and toxicities of these compounds as indicators. In addition, Wang et al. [15] quantitatively evaluated variations in the indoor air concentration of trichloroethylene (TCE) due to soil vapor intrusion into houses located both above and adjacent to the contamination plume in an aquifer. Their work indicated that houses that are laterally offset from the contamination plume are less affected by vapor intrusion because of the minimal transverse flux of TCE within the plume. However, the majority of such studies have considered only the steady-state flow behavior of the gas phase in the unsaturated zone, using constant concentration of VCSs as a boundary condition, and have not considered changes in the saturation distribution within the pore spaces. Consequently, quantitative predictions of unsteady variations in the volatilization flux near the ground surface based on the permeation of rainfall into the surface soil and changes in the contaminant distribution depending on time and space have not yet been reported.
In the present study, we developed the prediction model of volatilization flux from unsaturated soil contaminated by VCSs, based on mathematical modeling of a series of phenomena under unsaturated conditions. These factors included the multi-phase flow consisting of non-aqueous phase liquids (NAPLs), water, and gases in conjunction with the permeation of rainfall into the surface soil, the volatilization/condensation of VCSs, and the adsorption/desorption of VCSs. Using this model, we assessed the mechanism by which the gas volume flux (GVF) is generated at the ground surface due to variations in the ambient weather. In addition, the correlation between the transport phenomena of VCSs in surface soil and the mass flux were quantitatively evaluated and the effect of rainfall on the generation of volatilization flux was clarified. By adopting this new model, it is possible to predict the volatilization flux occurring due to transfer of VCSs from soil to the atmosphere considering natural phenomena such as rainfall at the contaminated site as a realistic phenomenon, comparing with a conventional conceptual model. This paper describes in detail the newly developed predictive model and the results obtained from this approach, based on a numerical analysis.

Governing Equation
A mathematical model was constructed based on the threephase flow of a gas/NAPL/water system resulting from the discharge of an undiluted solution of VCSs defined as NAPL and permeation of rainfall into the surface soil. The goal is to achieve a quantitative evaluation for the generation of the volatilization flux at the ground surface and the transfer of VCSs into the aquifer by considering volatilization/condensation, elution/precipitation and partitioning between pore water and soil particles due to adsorption/desorption. The relationship between each phase and each component of the VCSs in this mathematical model is summarized in Fig. 2. The mathematical model is  The relationship between  each phase and each component  in the present mathematical  model for soil and groundwater contamination by volatile  chemical substances (VCSs) constituted from three-phase including gas, NAPL, and water phases, and each component also considering the presence of water and air in addition to VCSs. The porous media flow of each phase is primarily determined by Darcy's law, and the change of composition in each phase is defined using mole fractions ( w g,k x w,k and y n,k ) as variables. The mass conservation equations for each phase and each component are defined as follows.
The equation for the NAPL phase (corresponding to an undiluted solution of the VCSs) is The equation for the water phase is The equation for the gas phase is The equation for the VCS components in the NAPL phase is The equation for the dissolved components in the water phase is The equation for the gaseous components in the gas phase is The relationship between each saturation defined as the volume ratio in the pore fluid, is expressed as In Eqs. (4) to (6), suffix k indicates component number and, as shown in Fig. 2, k values of 1, 2, and 3 corresponded to the VCSs, water, and air, respectively. The sum of the mole fraction in each phase ( w g,k x w,k and y n,k ) is always equal to 1. Because the NAPL phase corresponded to a one component system constituted only from VCSs and with k equal to 1, y n,1 is also always 1, meaning that ∇ ⋅ w g,k ⋅ Kk rg g g ∇Φ g + ∇ ⋅ D g,k ∇ w g,k g S g + R vnc,k + R vwc,k + R vsc,k = w g,k g S g t .
In addition, the dissolution of air into the water phase is not considered in the numerical analysis. Therefore, the dissolved components in the water phase comprised the VCSs and water, such that k is equal to 2 and The gas phase consists of three components (the VCSs, water and air), meaning that k equals 3 and In this mathematical model, the saturations values ( S g , S n and S w ), mole fractions ( w g,k x w,k and y g,k ) in each phase are initially estimated based on flow calculations. Following this, the elution/precipitation and adsorption/desorption of these components are considered in order to update each value in the same time step based on solubility and distribution coefficients. Therefore, the terms for elution and adsorption are not contained in these conservation equations. These treatments are described in the Sect. 2.4.

Formulation of Relative Permeability in a Three-Phase System
Relative permeability (RLP) is one of the most important parameters related to the evaluation of the unsaturated permeability characteristics under the co-existence of multi-phase in a porous medium. As discussed, under the unsaturated condition contaminated with VCSs, threephase flow consisting of an undiluted solution of VCSs (defined as the NAPL), gases containing volatilized VCSs, and water need to be considered for the numerical analysis. The RLP values for three-phase system have been extensively studied in the fields such as oil reservoir engineering and the evaluation of contamination by petroleum hydrocarbons, and several models for numerical analysis have been proposed [19][20][21]. All such models are based on the extension of the RLP values for two-phase NAPL-water systems ( k rn, nw S w , k rw, nw S w ) and gas-NAPL systems ( k rg, gn S g , k rn, gn S g ), because it is challenging to directly determine RLP values for three-phase systems. In such models, the relations of k rg, gn S g = k rg S g and k rw, nw S w = k rw S w are directly applied to analysis of three-phase systems. On the other hand, when k rl S g is defined as the RLP to the liquid phase consisting of both NAPL and water, the relation of k rn, gn S g =k rl S g is assumed. The RLP to the NAPL phase (8) y n,1 = 1 .
in a three-phase system, k rn S w , S g , can be obtained by extending the estimated values for a two-phase system. As an example, Aziz and Settari [22] normalized the model for estimating k rn S w , S g proposed by Stone [20] in the form of In the above equation, when k 0 rn and k 0 rl are defined as the RLP to NAPL at irreducible water saturation S wi in an NAPL-water system ( = k rn,nw S wi ) and the RLP for the liquid phase at residual gas saturation S gr ( = k rl S gr ), respectively, the assumption of k 0 rn = k 0 rl is also applied. Setting k rw S wi = 0 , k rg S gr = 0 and k rn S wi , S gr = k 0 rn at S w = S wi and S g = S gr , k rn S w , S g as in Eq. (11) permits the analysis of a three-phase system as a function of S g and S w . In addition, if the flow conditions in the pore spaces are assumed as single-or two-phase systems, k rn S w , S g can be categorized according to three scenarios, as follows.
1. If k 0 rn = 1 is assumed at S w = S wi and S g = S gr , Eq. (11) can be simplified as In the case of single-phase flow of a gas phase ( S l = 0 ), k rg (1) = 1 and k rl (1) = 0 must be satisfied. Substituting k rn,nw (0) = 1 and k rw (0) = 0 under the condition of S w = 0 corresponding to S l = 0 , k rn (0,1) = 0 can be estimated as 2. If there is no water phase in the pore spaces (that is, S w = 0 ), k rn, nw (0) = 1 and k rw, nw (0) = 0 must be satisfied. In this case, since k rn 0, S g for a three-phase system becomes identical to the RLP value for a liquid phase, we can obtain the relation of k rn 0, S g = k rl S g 3. In the case of S g = 0 corresponds to a two-phase system consisting of NAPL and water, substituting k rl (0) = 1 and k rg (0) = 0 into Eq. (12), k rn S w , 0 can be considered equivalent to k rn, nw S w and the relation of k rn S w , 0 = k rn,nw S w can be obtained as In the RLP model proposed by Stone [20] based on Eqs. (13) to (15), a consistent limitation for estimation of k rn S w , S g is applied even if phase condition in the pore spaces transitions from single-to three-phase. As noted, the RLP value for water, expressed as k rw S w , is assumed to be solely a function of S w even in three-phase system. Since k rw (0) automatically becomes 0 in a single-phase system comprising only a gas phase (meaning that S l = 0 ) and in a two-phase gas-NAPL system ( S w = 0 ), there is no contradiction in the treatment of water flow in the numerical analysis. In addition, the estimated value of k rw S w at S g = 0 can be directly applied to an NAPL-water twophase system. However, when the phase condition transitions from three-phase to a gas-water two-phase system, k rw S w does not change continuously because multi-phase flow conditions on the k rw S w curve are defined by S wi and by extent of the residual NAPL saturation S nr . In this case, as with the limitation for k rn S w , S g as expressed in Eq. (14), it is necessary to have k rw S w identical to k rl S g by treating as a function of both S g and S w . The flow behavior under the unsaturated condition that is the focus of the present study is primarily dominated by gas-water two phase flow with permeation of rainfall except ongoing contamination caused by discharge of an NAPL. Therefore, the treatment of permeability characteristics of water under three-phase condition is very important for a numerical analysis. In the present study, while maintaining the concept for the estimation of k rn S w , S g by Stone [20] shown in Eqs. (11) to (15), based on RLP curves for a gas-liquid two-phase flow, we attempted to formulate an RLP for a three-phase system by considering the continuity of RLP changes due to transitions in the phase condition. Figure 3a presents an example of RLP curves plotted as functions of S g for a gas-liquid two-phase system. In this mathematical model, k rg S g and k rl S g for the multiphase flow condition defined by S gr and the residual liquid saturation, S lr , are estimated based on Brooks and Corey's model [23], using the equations and (15) Here, gl is the pore size distribution index for the gas-liquid system and k 0 rg and k 0 rl are RLP values at each residual saturation ( S g = S gr and S g = 1 − S lr ), which are assigned arbitrary values less than 1. In Eqs. (16) and (17), S lr is assumed as the sum of S wi and S nr as Because Eqs. (16) and (17) are available only in the saturation interval defined by S gr ≤ S g ≤ 1 − S lr , corresponding to multi-phase flow, the applicable range of RLP curve must be extended by using k 0 rg and k 0 rl that correspond to S gr and S lr .
1. In the range of 0 ≤ S g ≤ S gr , the gas phase is retained inside the pore spaces and it always keeps k rg S g = 0 , while k rl S g will have a maximum value of k 0 rl under the multi-phase flow condition at S g = S gr . In addition, k rl S g = 1 at S g = 0 must be satisfied because this corresponds to the single-phase flow of liquid. Therefore, we applied the limitation that k rl S g changes linearly from 1 to k 0 rl as a function of S g in the range of 0 ≤ S g ≤ S gr , where k rg S g and k rl S g are defined as and (18) S lr = S wi + S nr .
In the range of 1 − S lr ≤ S g ≤ 1 , the liquid phase without mobility is retained inside the pore spaces because of S l < S lr and it always keeps k rl S g = 0 , while k rg S g will have a maximum value of k 0 rg in multi-phase flow condition at S g = 1 − S lr . In addition, k rg S g = 1 at S g = 1 must be satisfied. Therefore, we applied the limitation that k rg S g changes linearly from k 0 rg to 1 as a function of S g in the range of 1 − S lr ≤ S g ≤ 1 . k rg S g and k rl S g are defined as and As discussed, in order to maintain continuity of RLP curves along with the phase transitions from a three-phase to two-phase system, the RLP to water ( k rw S w , S g ) in a three-phase system must be treated as a function of both S g and S w as well as that for NAPL ( k rn S w , S g ). This is necessary because a multi-phase flow condition between NAPL and water will vary depending on the presence of (20) Formulation of relative permeability in a gas-NAPL-water three-phase system. a Relative permeability curves for a gas-liquid two-phase system. b Relative permeability curves for an NAPL-water two-phase system a gas phase in the pore spaces. In addition,k rn S w , S g at S w = 0 and k rw S w , S g at S w = 1 corresponding to twophase conditions consisting of gas-NAPL and gas-water must become identical to k rl S g , expressed as Eqs. (17), (20), and (22). In this case, as demonstrated in Fig. 3b, the value of k rl S g at a given S g value in a gas-liquid system must be directly introduced into the endpoint (with S w = 0 for NAPL and S w = 1 − S g for water) on the RLP curve for NAPL-water system. In addition, when k 0 rn and k 0 rw are respectively defined as the RLP values at each residual saturation ( S w = S wi and S w = 1 − S nr ) in an NAPL-water system, k rn S w , S g and k rw S w , S g at S w = S wi and S w = 1 − S nr − S g in a three-phase system can be obtained by multiplying k 0 rn and k 0 rw by k rl S g . As well as gas-liquid system as above, if Brooks and Corey's model [23] is applied to a multi-flow condition between NAPL and water, k rn S w , S g and k rw S w , S g will be estimated according to the following classifications depending on the change in the S w value in the pore spaces.
1. In the case of 0 ≤ S w ≤ S wi , while it always keeps k rw S w , S g = 0 , k rn S w , S g is estimated based on a linear equation including ( 0 , k rl S g ) and ( S wi , k rl S g ⋅ k 0 rn ), written as and 2. In the case of S wi < S w < 1 − S nr − S g , the estimations of k rn S w , S g and k rw S w , S g for multi-phase flow condition corresponding to the range from S w = S wi to S w = 1 − S nr − S g is directly based on Brooks and Corey's model [23] involving the equations and 3. In the case of and written as (27) k rn S w , S g = 0 Table 1 shows the classification system used to estimate k rn S w , S g and k rw S w , S g based on Eqs. (23) to (28). In this table, all the saturation combinations consisting of gas, NAPL and water phases supposed for the estimation of the RLP values are described. As an example, because k rl S g always has a value of 1 based on Eq. (20) in an NAPL-water twophase system at S g = 0 , the estimated values of k rn S w , S g and k rw S w , S g obtained from Eqs. (23) to (28) can be directly applied to NAPL-water two-phase system. In contrast, S g ≥ 1 − S lr , k rl S g is always set to 0 because the liquid phase does not have mobility. As a result, both k rn S w , S g and k rw S w , S g become always 0 for all saturation condition supposed. Therefore, the RLP model focusing on a three-phase system proposed in this study is able to consistently establish the continuity of RLP curves in association with transitions from three-phase to two-phase or from two-phase to single-phase.
As shown in Eqs. (16) to (28), because the residual saturations of each phase ( S wi , S nr , S gr and S lr ) define the interval of saturation corresponding to multi-phase flow condition on the RLP curves, these residual saturations must be quantitatively determined through experimental data prior to formulating the RLP. In addition, it is very important to consider the effect of particle size on the changes of these residual saturations, because the distribution of soil particle sizes varies widely. Based on S wi , S nr , and S gr values obtained from a series of soil column tests for NAPL-water and gas-water systems, we formulated as a function of the average diameter D A (m) of soil particle [24,25]. Figure 4 summarizes the effect of particle size on the changes in S wi , S nr , and S gr . Here, the horizontal axis shows the ratio of the particle diameter to the diameter for Toyoura sand, D A0 (= 2.00 × 10 −4 m), ( D A ∕D A0 ) which was used as a standard condition. From this figure, it was evident that the S wi value obtained from an NAPL-water system was highly dependent on the soil particle size. Specifically, as D A ∕D A0 decreased from 1.00 to 0.443, S wi increased from 0.186 to 0.491. This result is attributed to several effects. Firstly, smaller particles will have a larger specific surface area and thus have greater tendency of water-wet. Secondly, the effect of capillary pressure is increased in the case of smaller particles, such that more water is retained in the pores. On the other hand, the effect of particle size on S nr was not confirmed because the wettability of the pore spaces for NAPL was apparently intermediate between the wettability for gas and water. Even though the value of S nr is kept constant, increases in S wi decreases the value of the RLP to water. Therefore, it is interpreted that an NAPL flows more readily relatively with decreases in the particle size. Based on the experimental data, S wi can be estimated for the D A range from 0.00 to 2.00 × 10 −4 m based on the S wi0 (= 0.230) obtained for Toyoura sand as a standard value, using the equation Only water (25) Gas-water (27) Only gas S nr is assumed to be constant not affected by the soil particle size, and the average value of S nr0 (= 0.217) obtained for each condition is applied to the numerical analysis such that S gr showed a tendency to become lower with the decreases in D A ∕D A0 ratio, because gas phase behaves as non-wettable one in pore space. This result indicates that gas phase flows relatively easier with the decrease of particle size comparing with liquid phase. Based on the experimentally obtained values ranging from 0.039 to 0.287 for different D A , S gr is estimated using the S gr0 (= 0.287) for Toyoura sand as a function of D A as As shown in Eq. (18), the residual liquid saturation, S lr , in a gas-liquid two-phase system is assumed to equal the sum of S wi and S nr as expressed by Eqs. (29) and (30). Figure 5 provides a comparison of ternary diagrams summarizing the relative permeability values for the gas-NAPL-water three-phase systems with D A values set to (a) 2.00 × 10 −4 m and (b) 5.00 × 10 −5 m. The residual saturation ( S wi , S nr , S gr , and S lr ) determined from Eqs. (18) and (29) to (31) are also included in this figure. In addition, the pore size distribution index values ( gl , nw ) and RLP value at the residual saturations ( k 0 rg , k 0 rl , k 0 rn and k 0 rw ) are all set to 2.00 and 0.400, respectively. Each vertex in (30) S nr = S nr0 . Irreducible water saturation Residual NAPL saturation Residual gas saturation Fig. 4 The effect of soil particle size on variations in the residual saturation values S wi , S nr , and S gr Fig. 5 The ternary diagrams of relative permeability for gas-NAPL-water three-phase systems having average soil particle diameters, D A , of a 2.00 × 10 −4 m and b 5.00 × 10 −5 m the ternary diagram indicates a single-phase flow condition for each phase, and the value of RLP at the vertex is 1. The variations in the flow conditions as a function of the RLP for a three-phase system can be classified into eight zones, and the area labeled (3) in this figure represents the region within which the values of RLP for all phases are greater than 0. This corresponds to the range of saturation for a three-phase flow. This figure also indicates that, as D A becomes smaller, the three-phase flow condition shifts toward lower S n and higher S w values.

Treatment of Capillary Pressure in a Three-Phase System
In the present numerical analysis of a three-phase system consisting of gas, NAPL, and water, the system pressure, P (Pa), in the governing equations is generally considered to equal that of the gas phase, P g (Pa). The relationships between the pressure of each phase and the operating capillary pressure (CP) are and for the gas, NAPL, and water phases, respectively. Here, P n is the pressure of the NAPL phase (Pa), P w is the pressure of the water phase (Pa), P c,gn is the CP between the NAPL and gas phases (Pa) and P c,gw is the CP between the gas and water phases. If the CP between the NAPL and water phases is defined as P c,nw (Pa), the value of P c,gn in Eq. (33) can be obtained as The P c,gw value for a gas-water system is estimated using Brooks and Corey's model [23] based on the change in liquid saturation as In Eq. (36), S * l is the normalized liquid saturation and can be classified by considering hysteresis in the CP curve as or (32) P g = P, In the case that S g is set to 0 in Eq. (37), the CP curve is defined by the saturation interval between S l = S lr and S l = 1 . This condition corresponds to the first drainage process by gas as a non-wettable phase into the pore spaces saturated with water. On the other hand, the estimation of P c,gw follows Eq. (38) in the case of S l < 1 − S gr , the obtained CP curve corresponds to the imbibition process of water as a wettable phase. Whereas variations in S g during both the drainage and imbibition process are considered for changes in the CP curve in the case of 1 − S gr ≤ S l ≤ 1 , the curve obtained for S l < 1 − S gr is treated as constant shape even if S g varied. The value of P d,gw in Eq. (36) is defined as the threshold pressure (Pa) when the permeation of gas into the pore spaces saturated with water is initiated, and this value can be estimated using the Young-Laplace equation Here, gw is the surface tension in a gas-water system (N/m), is the contact angle of water (rad.) and d pore is the principal pore diameter (m). When the hysteresis of the CP curve is considered, the limitation for which P c,gw becomes 0 at a specific S g with S l = 1 − S g + S glimit as the threshold can be applied in the case of S l ≥ 1 − S gr (as shown in Fig. 6a) based on the equation Here, N gw is reduction order of P c,gw . Figure 6b shows a comparison of CP curves at 20 °C for different D A , in which the solid and dashed lines indicate CP curves for the first drainage and imbibition processes, respectively. In this figure, the residual saturation ( S gr , S lr ) in Eqs. (37), (38), and (40) were estimated by using Eqs. (18) and (29) to (31) and lg , S glimit and N gw were set to 2.00, 0.05 and 4.00, respectively. In addition, gw and in Eq. (39) were given values of 7.27 × 10 −2 N/m and 1.23 rad based on literature sources [26,27]. The values of d pore are calculated using D A value obtained from a nuclear magnetic resonance analysis of a core sample via the equation [28] As shown in Fig. 6b, as D A becomes smaller, the P c,gw value for a given S l condition increases, due to the higher S lr and smaller d pore , based on relationship Here, P d,nw is threshold pressure at the initiation of NAPL permeation as a non-wettable phase and S * w is the normalized water saturation. S * w can be classified into two cases in order to maintain the continuity of the CP curve with the transition from three-phase to two-phase conditions and considering the hysteresis of the CP curve in the case that S n ≤ S nr . The associated equations are and Furthermore, in the case that S w ≥ 1 − S nr − S g , introducing S nlimit as the threshold, the limitation for which P c,nw becomes 0 at S w = 1 − S n − S g is applied for the estimation of P c,nw via the equation Here, N nw is the reduction order of P c,nw . Based on Eqs.
(43) to (45), while the change in P c,nw due to both S g and S n is considered if 1 − S n − S g ≤ S w ≤ 1 − S g , P c,nw changes depending only on S g if S w < 1 − S nr − S g .

Treatment of the Volatilization and Condensation of VCSs
The quantitative evaluation of the volatilization flux at the ground surface required the effect of the volatilization/ (43) condensation of VCSs to be considered in conjunction with calculations of flow in the surface soil. In this mathematical model, the total volatilization (or condensation) rate for each phase, expressed as R vn , R vw and R vs (kmol/m 3 /s) as in Eqs. (1) to (3), is defined as the sum of each component ( R vnc,k , R vwc,k and R vsc,k ) and introduced into the mass conservation equations expressed as Eqs. (1) to (6). The volatilization rate for the NAPL phase/condensation rate from gas phase is determined as while the volatilization rate from the water phase/condensation rate from the gas phase is and the volatilization rate derived from the adsorbed VCSs on the soil particle surfaces is The volatilization of each component from the liquid phase and partitioning between the gas and liquid phases in the pore spaces are generally dependent on the saturated vapor pressure, P sat,k (Pa). Based on the Clausius-Clapeyron equation, P sat,k at a given temperature, T (K), can be determined from the saturated vapor pressure, P sat0,k (Pa), at a standard temperature, T s (K), and the heat of vaporization, ΔH v,k (J/kmol), as In addition, the mole fraction, y l,k (-), of component k in the liquid phase containing VCSs adsorbed on soil particles is expressed as Here, x s,k is the adsorption concentration of component k (kmol/kg). According to Raoult and Dalton's law, the saturated vapor pressure of each component in a multi-component system can be expressed as the product of y l,k and P sat,k , and equals the partial pressure, w g,k P , of each component in the gas phase at vapor-liquid equilibrium. Assuming that each component achieves vapor-liquid equilibrium in the pore spaces ( y l,k P sat,k = w g,k P ), the rate of volatilization (or condensation) per mole of component k, R v0,k (1/s), can be expressed as Here, k v is the rate constant necessary to reach gas-liquid equilibrium [1/(Pa s)]. The volatilization/condensation rates ( R vnc,k , R vwc,k and R vsc,k ) of component k in each phase can be classified as either y l,k P sat,k > w g,k P or y l,k P sat,k < w g,k P . In the case of y l,k P sat,k > w g,k P , component k volatilizes from the liquid phase and is distributed into the gas phase, such that and On the other hand, in the case of y l,k P sat,k < w g,k P , component k condenses from the gas phase and is distributed into the liquid phase. The partition between the NAPL and water phases with condensation depends on the molar ratio of component k in the liquid phase as and

VCS Elution into the Water Phase and Adsorption to Soil Particles
For the prediction of distribution of VCSs in the surface soil depending on time and space, the elution of these compounds into the water phase and adsorption onto soil particles are very important parameters. In this mathematical model, it (50) y l,k = (yn,k n S n +x w,k w S w) +x s,k s (1− ) n S n + w S w .
(51) R v0,k = k v y l,k P sat,k − w g,k P .
(52) R vnc,k = R v0,k y n,k n S n , (55) R vnc,k = y n,k n S n y n,k n S n + x w,k w S w R v0,k w g,k g S g , is assumed that the equilibriums of elution and adsorption between NAPL, pore water and soil particles are simultaneously established and that the time step for the calculation of fluid flow is sufficiently long, comparing with the time required to reach the equilibrium conditions for elution and adsorption. Based on these assumptions, the elution of the VCSs into the water phase and adsorption onto the soil particles are treated as follows. Initially, the mole fractions for the NAPL and water phases ( x w,k and y n,k ) at each location are calculated by solving advection-dispersion equations [see Eqs. (4) and (5)], after which x w,k , y n,k and the adsorption concentration, x s,k , are updated in the same time step based on the octanol/water partition coefficient, K ow,k (m 3 /m 3 ), and the distribution coefficient between soil and water, K d,k (m 3 /kg).
Expressing the block lengths in each direction as Δx , Δy and Δz (m), the total moles, N mol,k (kmol), of each VCS component k per block are calculated as Since the mole fraction ( y n,k , x w,k ) is converted into a concentration in each phase (kmol/m 3 ) through multiplying by the molar density ( n , w ), the elution equilibrium between the NAPL and water phases can be obtained using K ow,k as In addition, K d,k and x w,k are related to x s,k according to the equation By substituting the right-hand sides of Eqs. (59) and (60) into Eq. (58), the value of x w,k can be updated taking into account the effects of dissolution and adsorption as If the x w,k value estimated from Eq. (61) exceeds the solubility value x w,sat,k , x w,k is updated to x w,sat,k . Using the resulting x w,k , y n,k and x s,k are also updated based on Eqs. (59) and (60). If x s,k exceeds the amount of saturated adsorption, x s,sat,k , x s,k is again updated, to x ssat,k and the residue is re-distributed into the NAPL phase.

Validation of Numerical Analysis Model
To demonstrate the validity of this mathematical model, we had conducted history matching for the occurrence behavior of volatilization flux observed in the laboratory-scale experiment using soil column [29] (Kondo et al., submitted for publication). Diameter and length of soil column were 7.0 cm and 27 cm, and porosity was 0.42. Soil column was set up inside (58) N mol,k = y n,k n S n + x w,k w S w + s (1 − )x s,k ΔxΔyΔz.
(59) n y n,k = K ow,k w x w,k .
(60) x s,k = K d,k w x w,k .
water bath vertically, and temperature of soil column was controlled at a given value by circulation of refrigerant. In this experiment, ethanol (C 2 H 5 OH) was used as model substance of VCSs. As the experimental procedure, firstly, ethanol solution of 10,000 mg/L was injected from the bottom edge of soil column under dry condition. The injection of ethanol solution was terminated when the increase of ethanol concentration due to the arrival of the front of injected ethanol at the top edge of column was detected by the gas sensor installed at the outlet of column. After that, based on the concentration change of ethanol in gas phase at the outlet of column by changing column temperature from 42 to 12 °C, the occurrence behavior of volatilization flux of ethanol depending porous media flow under unsaturated condition, advection and dispersion in gas and water phases was analyzed, and obtained series of data for the temperature dependence on the occurrence of volatilization flux. Figure 7 shows schematic of analytical mesh zone and boundary condition corresponding to the apparatus of soil column test (left) and the comparison of experimental and calculated results for the changes of temperature and volatilization flux of ethanol with time (right). For history matching process, we changed relative permeability to gas and water phase ( k rg , k rw ), volatilization rate constant k v and dispersion coefficient D g,k as parameters. Through the optimization of these parameters, the calculated result for the changes of temperature and volatilization flux of ethanol could reproduce the experimental ones sufficiently, we have confirmed the validity of this mathematical model.

Analytical Mesh Zone and Boundary Conditions
In this mathematical model, the vertical transfer of fluid and contaminants resulting from the permeation of rainfall and the effect of capillary pressure was dominant in the surface soil, and so each mass conservation equation [Eqs. (1) to (6)] was discretized in the x-z two-dimensional coordinate system. As a solution for the numerical analysis, after discretization of each conservation equations by using the finite difference method, implicit method for pressure and an explicit method for saturation and concentration were applied as the solving methods, respectively.
The x-z two-dimensional analytical mesh zone employed in this study is presented in Fig. 8. Boundary blocks corresponding to air on the ground were arranged at the top edge of the mesh zone. The porosity, , and absolute permeability, K , in the horizontal direction for these blocks were set to 1.00 and 98.7 μm 2 (= 100 Darcy), respectively. Each residual saturation ( S wi , S nr , S gr , and S lr ) was given a value of 0 and the effect of capillary pressure was not considered so as to facilitate fluid flow compared to the inside of the pore spaces. The changes in temperature, air pressure and rainfall level as summarized in Fig. 8 were applied to these blocks as boundary conditions and only gas phase diffusion was considered in the case of the boundary surface at the top edge. In addition, it was assumed that gaseous components derived from volatilization present in these blocks diffused instantaneously into the atmosphere. Therefore, the mole fractions of the VCSs and of the water vapor ( w g,1 , w g,2 ) were set to 0 for the estimation of the rate of volatilization (or condensation) per mole of component k, R v0,k [1/s] using Eq. (51). Each block at the bottom edge of the mesh corresponded to the boundary with an aquifer. The water saturation, S w , at the boundary surface with the aquifer was set to a value of 1 and a change in air pressure equal to the ground level was applied to each block as a boundary condition, since the first aquifer located just below the surface soil was not pressurized. Seasonal variations in the temperature distribution within the surface soil were taken into account based on the temperature gradient between the ground surface and the water table, while the temperature of the aquifer was held constant at 15 °C. The lateral boundary conditions were decided through preliminary calculations consisting of the three steps shown in Fig. 9, and these conditions were subsequently introduced into the numerical analysis. In step 1, the steady-state distributions of pressure and saturation considering only rising of the water table due to suction without rainfall were estimated. In Step 2, the distributions of pressure and saturation obtained in step 1 were applied as constant lateral boundary conditions and calculations were carried out by considering weather changes as a boundary condition for each block corresponding to air on the ground. Finally, in step 3, the seasonal variations in the Fig. 8 The analytical mesh zone and boundary conditions for the x-z two-dimensional system employed in the numerical analysis vertical distributions of pressure and saturation at the central part obtained in step 2 were introduced as the lateral boundary conditions during the actual computational process. Table 2 summarizes the basic parameters of surface soil and the physical parameters for benzene as a model VCS contaminant in the numerical analysis. Based on the analytical mesh zone shown in Fig. 8, the lengths in x and y directions and the thickness in the z direction were set to 10.0, 1.00, and 4.00 m, respectively. The weather conditions applied as boundary conditions to the air on the ground, including temperature, air pressure and rainfall level, were based on a series of data acquired in Ibaraki prefecture in 2017 as reported by the Japan Meteorological Agency [30], and introduced for 10 years calculation in one year cycles. D A was set to 5.00 × 10 −5 m, corresponding to silt, and the absolute permeability, K , was estimated to be 2.95 × 10 -2 μm 2 based on the work of Matsuo and Kogure, who reported the equation [31] Here, w,15 • C is viscosity of water at 15 °C (= 1.14 × 10 −3 Pa s). The initial total petroleum hydrocarbon (TPH) was set to
1.00 × 10 4 mg/kg and the depth range initially contaminated by benzene was set from − 0.10 to − 1.00 m from ground level. According to the relationship between each phase and each component in this mathematical model (shown in Fig. 2), the 3.53 × 10 2 component number, k, for benzene in the multi-component system was defined as 1. As physical parameters related to elution and adsorption, the octanol/water partition coefficient, K ow,1 , the distribution coefficient between soil and water, K d,1 , and the solubility of benzene into water phase, C sat,1 , were set to 1.35 × 10 2 m 3 /m 3 , 5.89 × 10 −2 m 3 /kg, and 1.75 × 10 3 mg/L at 25 °C, respectively [32]. In addition, the saturated vapor pressure, P sat0,1 , at 20 °C and standard boiling point, T b,1 , were set to 1.00 × 10 4 Pa and 3.53 × 10 2 K based on values in the literature [33]. Below, we describe the calculation results and discuss the mechanism associated with the generation of GVF at the ground surface as well as the effect of each parameter on the magnitude of the benzene mass flux due to volatilization.

Mechanism for the Generation of the GVF at the Ground Surface
Here, we initially consider the mechanism by which GVF is generated at the ground surface due to weather variations. The GVF is defined as the total volume flow rate (SL/ m 2 /day) for the gas phase containing the volatilized VCSs, water vapor and air per unit area of the ground surface.
This parameter can be estimated based on Darcy's law according to the upward pressure gradient of the gas phase at the ground surface. Even if the volatilization of a VCS component such as benzene occurs in the pore spaces, it is supposed that this volatilized component does not diffuse from the ground surface to the air so much due to the effect of advection since the amount of generated GVF becomes relatively small in the case that the gas permeability is relatively low. Therefore, the magnitude of the GVF is a very important parameter related to the quantitative evaluation of the volatilization flux of a VCS component. Figure 10 shows the effect of weather changes on the GVF based on calculation result obtained for the first year of a 10-year time span. It is apparent that the GVF increased from summer to autumn with temperature increase and the air pressure decrease. The peak GVF value in this case reached 1.58 SL/ m 2 /day at the 0.6 year point, and the cumulative GVF at the end of the 1-year period was estimated to be 0.206 Sm 3 /m 2 . These data reflect the promotion of volatilization of VCS and water due to increase of temperature and decrease of air pressure as well as increases in the upward pressure gradient at the ground surface due to due to the rise of pore pressure. These results exhibit good agreement with previous studies that attempted to clarify correlations between GVF values and temperature increases derived from solar radiation [34], and demonstrated dramatic increases in the GVF in conjunction with rainfall events [35]. In addition, as can be observed after the 0.8 year point, the presented data also show the characteristic behavior that the GVF significantly decayed and was maintained at a lower level for a certain period after a heavy rainfall. Thus, this model was able to reproduce the dynamic behaviors of VCSs, water vapor and air in soil at the boundary of the subsurface soil, such that the simulated trends were in good agreement with the monitoring data.
The changes in pressure gradient (atm), water saturation, S w (-), and relative permeability to the gas phase, k rg (-), at the ground surface ranging from 0.2 to 0.4 years in addition to the GVF are expanded in Fig. 11, to allow a discussion of the mechanism of GVF generation with the permeation of rainfall into the surface soil. These data indicate that the gas phase pressure in the pore spaces maintained higher than the ambient atmospheric pressure at the ground surface because volatilization was promoted so as to maintain saturated vapor pressure of each component existing in liquid phase. As such, pressure gradient always indicated positive value toward air on the ground except the period just after heavy rainfall. Therefore, it can be interpreted that the gas phase generally flows upward from the interior of the surface soil toward above the ground. In the period labeled (1) spanning from 0.27 to 0.30 years, a heavy rainfall with a total of 89.5 mm precipitation occurred and, as a result, S w increased from 0.936 to 0.972 due to the permeation of rainfall into the surface soil. During this time span, the GVF exhibited a rapid decay from 5.49 × 10 −1 to 0.00 SL/m 2 /day due to the decrease in k rg , since the pore spaces approached residual gas saturation S gr . In contrast, during period (2) from 0.30 to 0.35 years just after this heavy rainfall, S w gradually decreased Fig. 11 The mechanism by which the gas volume flux is generated at the ground surface from 0.972 to 0.884, and it was confirmed that pore space had restored un-saturated condition. As a result, the increase in k rg due to the decrease in S w induced an increase in the GVF from 0.00 to 1.12 SL/m 2 /day. From these trends, it is apparent that the increase in k rg was the dominant factor affecting the generation and magnitude of a GVF at the ground surface. Figure 12 summarizes the changes in the distributions of TPH (mg/kg), NAPL saturation, S n (-), the concentration of benzene dissolved in the water phase, C w (mg/L), and the amount of benzene adsorbed on the soil particles, C s (mg/kg), in the vertical direction over time. Whereas the initial TPH concentration was set to 10,000 mg/kg equivalent to an undiluted solution of benzene defined as the NAPL phase, the subsequent distribution of TPH shown in Fig. 12 corresponds to the total amount of benzene at each position based on elution into the water phase and adsorption onto the soil particles, in addition to the undiluted solution.

Correlation Between Transport Phenomena of VCS Component in Surface Soil and Variations in the Mass Flux at the Ground Surface Due to Volatilization
Although the S n value corresponding to the initial TPH was 0.04, the amount present in the pore spaces as the NAPL phase after taking into account partitioning between the water and solid phases based on Eq. (58) was estimated to be approximately 0.015 in terms of S n . In this case, since pore space was under the condition below residual NAPL saturation S nr based on Eq. (30), benzene as the undiluted solution did not have mobility depending on relative permeability to NAPL phase k rn . From this figure, it is apparent that S n became lower gradually beginning from the ground surface over time as elution of benzene into water phase proceeded due to the permeation of rainfall. These data confirm that the elution process was complete after 5 years when S n had reached 0. In addition, whereas a high C w value was maintained (corresponding to a solubility of 1.75 × 10 3 mg/L) as long as the undiluted benzene remained in the pore spaces, C w after the completion of the elution process decreased to approximately 1.30 × 10 2 mg/L and benzene in the water phase spread downwards due to the effect of advection and dispersion while maintaining this concentration level. The saturated adsorption of benzene, x s,sat,1 , in this case was tentatively set to 1.00 × 10 −4 kmol/kg, corresponding to 7.81 × 10 3 mg/kg, and C s was maintained at this value regardless of variations in C w because benzene had high adsorptive onto the soil particles. The maximum TPH was maintained at 10,000 mg/kg, which was equal to the initial condition, until 2 years later when an undiluted solution of benzene remained present. Subsequently, the maximum TPH was estimated to be approximately 6000 mg/kg and the data show that the majority of the TPH was derived from adsorption because C w decreased. Based on the distribution of the TPH, we confirmed that the spread of benzene as a pollutant after 10 years did not exceed a depth of 2 m, reflecting the high adsorptive of benzene on the soil. Although the mobility of benzene in the surface soil was relatively low, as shown in Fig. 12, TPH at the ground surface decreased due to the progress of elution into the water phase over time. Figure 13 shows the changes in the TPH, GVF, mass flux of benzene (MFB) and cumulative MFB at the ground surface with time. Here, MFB is defined as the mass flow rate (mg/m 2 /day) of benzene in the gas phase per unit area of ground surface. This value was calculated as the sum of the benzene outflow resulting from both advection and dispersion from the ground surface into the air on the ground. The MFB value due to advection could be obtained by multiplying the GVF (m 3 /m 2 /day) by the molar density, g (kmol/m 3 ), the mole fraction of benzene in the gas phase, w g,1 (-), and the molar mass, M g,1 (kg/kmol), of benzene. The MFB value due to dispersion was based on the concentration gradient of benzene in the gas phase at the ground surface corresponding to ΔM g,1 w g,1 g S g (kg/m 3 ). From this figure, it is evident that the TPH at the ground surface decreased to 1800 mg/kg after 10 years from an initial value of 10,000 mg/kg. Since this decrease reduced the mole fraction of benzene, y l,1 , in the liquid phase (as defined by Eq. (50)), the saturated vapor pressure of benzene, P sat,1 , in the multicomponent system was also lowered so as to decrease the volatilization rate, R v0,1 , calculated as in Eq. (51). Consequently, whereas the GVF changed periodically depending Fig. 14 The correlation between transport phenomena of benzene in surface soil and the generation of mass flux of benzene (MFB) in association with a variations in temperature, air pressure and soil moisture content and b total petroleum hydrocarbon (TPH) after 1, 2, and 5 years only on weather variations, these data confirmed that the MFB decayed over time after MFB showed a maximum value of 114 mg/m 2 /day at 0.56 years. The cumulative MFB after 10 years was estimated to be 109 g/m 2 , equal to just 6.81 × 10 −1 % of the initial TPH. This result indicates that almost none of the benzene initially present was discharged from the ground surface as a volatilization flux. Rather, the most benzene remained in the pore spaces.
As discussed, increases in k rg as the pore spaces were restored to their unsaturated state due to decreased rainfall was the primary factor that determined both the generation and the magnitude of the GVF. Therefore, the diffusion of benzene into the air on the ground from the ground surface also depended significantly on the effect of weather changes on the transport phenomena in the surface soil. Figure 14 summarizes the correlation between the transport phenomena of benzene in the surface soil and the MFB generation, based on the temperature-air pressure-soil moisture content and temperature-air pressure-TPH concentrations interactions. From this figure, in addition to the effect of the increase in GVF due to temperature increase and air pressure decrease from summer to autumn, it was confirmed that the larger MFB was generated reflecting the increase of total amount of benzene existing in gas phase since gas saturation S g became higher in the case of lower soil moisture content. On the other hand, because the mole fraction of benzene in the gas phase, w g,1 (-), became higher when the TPH at the ground surface was relatively high during the initial stage of contamination, the MFB increased in the case of higher temperatures, lower air pressures and higher TPH. In particular, the results showed that temperature variations had the most significant effect on the magnitude of the MFB and that this value became much larger under the conditions of lower soil moisture contents and higher levels of TPH in summer season with high temperature.

The Effect of Rainfall Level on the Generation of Volatilization Flux at the Ground Surface
The effects of rainfall on the generation and magnitude of the volatilization flux were assessed by performing additional calculations in which the rainfall was set to 0, 600.5, or 2402 mm/year, and comparing the results to the value for a rainfall level of 1201 mm/year as a standard condition. In the case of the 600.5 and 2402 mm/year calculations, the amounts of rainfall shown in Fig. 10 were halved or doubled, respectively, and applied to the blocks corresponding to the air on the ground as a boundary condition. Figure 15 provides a comparison of the variations in the GVF and water saturation, S w , at the ground surface over the first year for different rainfall levels. In the case without rainfall, S w was maintained at lower level, ranging from 0.750 to 0.800. As a result, because the relative permeability to the gas phase, k rg , was increased, the GVF was more than five  16 The total petroleum hydrocarbon (TPH) and the cumulative mass flux of benzene (MFB) at the ground surface in the case of different rainfall levels times greater than that with rainfall, with a maximum GVF of 9.42 SL/m 2 ·day. On the other hand, as the rainfall level increased, S w was generally maintained at a higher level. The change of S w for each rainfall level ranged from 0.832 to 0.978 (1201 mm/year), from 0.822 to 0.972 (600.5 mm/year), and from 0.844 to 0.981 (2402 mm/year). In these cases that the average diameter of the soil particles, D A , was set to 5.00 × 10 −5 m, the residual gas saturation, S gr , was estimated to be 1.79 × 10 −2 based on Eq. (31). The decrease in k rg for a rainfall level of 2402 mm/year was especially remarkable since the period which S w showed close to a value corresponding to S gr became longer. Therefore, the GVF decreased as the rainfall level increased, giving maximum values of 1.53 SL/m 2 /day (at 600.5 mm/year), 1.38 SL/m 2 /day (at 1201 mm/year), and 1.12 SL/m 2 /day (at 2402 mm/year).
As shown in Fig. 16, the change in the MFB at the ground surface was correlated with the GVF and with variations in water saturation. This correlation can be derived as follows. In the case without rainfall, the effect of dilution on the dissolved component in the water phase was relatively small and the transportation of benzene in the surface soil was determined only by diffusion resulting from a concentration gradient, because the elution into the water phase and migration of benzene due to the permeation of rainfall did not occur. Consequently, the TPH at the ground surface gradually decreased compared to the other cases obtained with rainfall, and the TPH after 10 years was relatively high at 6500 mg/ kg. Because the GVF was increased and the mole fraction of benzene in the gas phase, w g,1 (-), was maintained at a higher level, there was a remarkable increase in the MFB. The cumulative MFB after 10 years was estimated to be 720 g/m 2 , which was several times the value obtained with rainfall. This result demonstrated that the concentration of VCSs in indoor air would be expected to increase if there is a lack of rainfall, leading to possible health risks related to the inhalation of these compounds. However, because the GVF decreased as the rainfall level increased, reflecting the decrease in k rg due to the increase in S w , and the total amount of benzene in the gas phase at a lower S g also became smaller, the cumulative MFB values were estimated to be 184 g/m 2 (at 600.5 mm/year), 109 g/m 2 (at 1201 mm/year), and 65.9 g/ m 2 (at 2402 mm/year). Therefore, the mass flux of VCSs at the ground surface generally decreased with greater rainfall, showing that the risk level in the open air would be relatively low due to the effect of diffusion in the atmosphere.

Conclusions
This study developed a model for predicting the volatilization flux of VCSs at the ground surface. Using this developed model, we consider the mechanism of generation of GVF at the ground surface due to weather variations. Targeting benzene as a model VCS, the correlation between transport phenomena of benzene in surface soil and the change in the mass flux of benzene generated was quantitatively evaluated. The effect of rainfall on the generation of volatilization flux was also clarified. The main conclusions of this work can be summarized as follows.
Based on mathematical modeling for multi-phase flow consisting of an NAPL, water, and gas under the unsaturated condition, we developed the prediction model of volatilization flux derived from VCSs at the ground surface by considering volatilization/condensation and adsorption/ desorption of VCSs component. Relative permeability and capillary pressure in a three-phase system were newly modeled and formulated by considering the continuity and hysteresis associated with phase transitions from three-phase to two-phase or from two-phase to single-phase system. In addition, factors related to weather, including rainfall level, temperature and air pressure, were introduced as boundary conditions in order to evaluate the effects of such factors on the generation and magnitude of the volatilization flux.
Simulation results demonstrated significant increases in the GVF as the pore spaces restored an unsaturated condition after a heavy rainfall. In addition, increases in the relative permeability to the gas phase, k rg , was the dominant factor determining the generation and magnitude of the GVF at the ground surface. The simulated results also show that (1) whereas the GVF changed periodically depending only on weather variations, the MFB decayed over time since the TPH decreased with permeation of rainfall into the surface soil. (2) This mass flux became much larger at lower levels of soil moisture content in addition to the effect of increase in calculated GVF at higher temperature and lower air pressure.
(3) In the absence of rainfall, increases in the GVF and high mole fractions of benzene in the gas phase produced a very large increase in the MFB. These results indicate that the concentrations of VCSs in indoor air can undergo a relative increase comparing with that in the open air with rainfall, such that there may be increased health risks derived from the exposure through the inhalation of these compounds.
Author Contribution All authors contributed equally to this article.

Conflict of Interest The authors declare no competing interests.
Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http:// creat iveco mmons. org/ licen ses/ by/4. 0/.