Influence of groundwater table fluctuation on the non-equilibrium transport of volatile organic contaminants in the vadose zone

Journal Pre-proofs Research papers Influence of groundwater table fluctuation on the non-equilibrium transport of volatile organic contaminants in the vadose zone Shengqi Qi, Jian Luo, David O'Connor, Xiaoyuan Cao, Deyi Hou PII: DOI: Reference:

S0022-1694(19)31088-1 https://doi.org/10.1016/j.jhydrol.2019.124353 HYDROL 124353

To appear in:

Journal of Hydrology

Received Date: Revised Date: Accepted Date:

24 July 2019 17 October 2019 11 November 2019

Please cite this article as: Qi, S., Luo, J., O'Connor, D., Cao, X., Hou, D., Influence of groundwater table fluctuation on the non-equilibrium transport of volatile organic contaminants in the vadose zone, Journal of Hydrology (2019), doi: https://doi.org/10.1016/j.jhydrol.2019.124353

This is a PDF file of an article that has undergone enhancements after acceptance, such as the addition of a cover page and metadata, and formatting for readability, but it is not yet the definitive version of record. This version will undergo additional copyediting, typesetting and review before it is published in its final form, but we are providing this version to give early visibility of the article. Please note that, during the production process, errors may be discovered which could affect the content, and all legal disclaimers that apply to the journal pertain.

© 2019 Published by Elsevier B.V.

Influence of groundwater table fluctuation on the non-equilibrium transport of volatile organic contaminants in the vadose zone Shengqi Qi a, Jian Luo b, David O’Connor a, Xiaoyuan Cao a, Deyi Hou a,*

a

School of Environment, Tsinghua University, Beijing 100084, China

b

School of Civil and Environmental Engineering, Georgia Institute of Technology, Atlanta, GA 30332-0355, United States * Corresponding author: Deyi Hou ([email protected])

Abstract The migration of volatile organic compounds (VOCs) from groundwater is an environmental concern because of the risk of vapor intrusion. In coastal aquifers or areas under the influence of systematic pumping, groundwater table fluctuations may be significant, which may affect VOC flux. This study developed a one-dimensional numerical model to simulate VOC mass transfer from groundwater to the vadose zone to the atmosphere, which was validated with data gained from a previously performed soil tank experiment. The numerical model revealed that groundwater table fluctuations in a simulated sandy aquifer cause VOC flux levels to increase by more than one order of magnitude. Furthermore, groundwater table fluctuations of larger amplitude are associated with larger VOC flux levels. Contaminants associated with larger Henry’s law constants and smaller water-solid distribution coefficients are influenced to a greater extent by groundwater table fluctuations. In a coarse sand aquifer, the nonequilibrium model demonstrated that the VOC flux level is lower (18%) than that calculated if equilibrium is assumed, showing the equilibrium assumption to be somewhat conservative. Groundwater table fluctuations induce significant increases in VOC flux from soils of coarse or medium sand, while flux levels are not significantly affected in silt or clay soils. Rainfall only has a temporal effect on VOC flux during groundwater table fluctuations. The numerical model in this paper could help improve estimates of human health risk by VI in areas with groundwater table fluctuation. Further non-equilibrium studies are needed to estimate mass transfer coefficients

among the three phases for different contaminants under different conditions.

Keywords: groundwater table; fluctuation; volatile organic compound; released flux; non-equilibrium model

1. Introduction Volatile organic compounds (VOC) in groundwater systems are a great environmental concern, with a growing imperative to deal with the hazard of widespread VOC plumes caused by anthropogenic pollution (Hou and Li, 2017). For instance, groundwater contamination involving the highly toxic VOC carbon tetrachloride (CCl4) (Recknagel, 1967) is a particular concern in China, due to its wide use by industry. The volatilization of VOCs from groundwater to the vadose zone represents an important migration pathway that threatens the health of humans in buildings located above (Garg et al., 2017; Rivett et al., 2011; Yao et al., 2018). Therefore, there is an essential need to formulate vapor intrusion models that can accurately predict VOC transport, so that health risks can be adequately managed. Many models have already been developed for this purpose, most of which conceptualize transport in the vadose zone by diffusion, advection, and biodegradation, assuming a static groundwater table and equilibrium partitioning in vapor, liquid and solid phases (DeVaull, 2007; Verginelli and Baciocchi, 2011). A retardation factor (Rt) is typically defined to integrate the mass distribution between gas phase and liquid phase, and between liquid phase and solid phase under equilibrium conditions (Costanza-Robinson et al., 2013; Kim et al., 2005; Shen et al., 2014; Shen and Suuberg, 2016). However, some scenarios may influence VOC transport that are not adequately addressed by these models (Cho et al., 1993). For example, time-varying groundwater systems involving different mechanisms such as pressure variations, rainfall events and

groundwater table fluctuations (Shen et al., 2012a; Shen et al., 2013c). Coastal aquifers around the world are subject to groundwater table fluctuations induced by oceanic oscillations (Shoushtari et al., 2016). Although the fluctuation amplitude would damp out with the distance from the shoreline (Raubenheimer et al., 1999), Singaraja et al. (2018) observed that the amplitude can be as much as 0.05~0.10 m within a distance of 700 m from the shoreline. In addition, Marechal et al. (2002) reported that groundwater tables may fluctuate with an amplitude of 0.03~0.10 m owing to tidal effects in some inland areas. Such fluctuation, especially the drainage process, facilitates the release of contaminants from groundwater to the vadose zone (Guo, 2015; Werner and Hohener, 2002). This factor could have profound implications for the world’s population living in close vicinity to VOC groundwater plumes in coastal areas. For instance, in the San Francisco Bay area of California, USA, there are >8,000 contaminated

sites

within

5

km

of

the

Bay

shoreline

alone

(http://geotracker.waterboards.ca.gov/). Although not all of these sites will involve VOC groundwater plumes, many of them will. Therefore, the influence of such scenarios in vapor intrusion models, and the validity of local equilibrium assumptions, needs to be further examined, with models considering non-equilibrium mass transfer possibly being a more realistic way to describe the transport. Some previous studies have considered the transport of VOCs from groundwater to the atmosphere under the condition of groundwater table fluctuations. However, they performed experiments with one or two cycles of groundwater table fluctuation

(Guo, 2015; Werner and Hohener, 2002). Also, some studies only considered the “block” effect of capillary fringe on the transport of VOCs from saturated zone to the unsaturated zone, while they did not specifically study the effect of groundwater table fluctuation (Shen et al., 2013a). Furthermore, most transport models assume equilibrium among different phases (Picone et al., 2012), while other models ignore the sorption of contaminants to the soil surface (Thomson et al., 1997). If steady state without groundwater table fluctuation is assumed, VOC transfer between gas and liquid phases, and between liquid and solid phases will be in equilibrium, and non-equilibrium vapor intrusion models with kinetic mass transfer would be simplified to equilibrium vapor intrusion models. As mentioned by Chen et al. (2016), mass transfer limitation could significantly change the distribution of gaseous concentration and liquid concentration during the process of vapor intrusion. However, there is a research gap regarding mass transfer limitation during groundwater table fluctuation. Therefore, the development of a comprehensive numerical model that considers non-equilibrium mass transfer during groundwater table fluctuations is needed. Until now, a number of models have been introduced to consider the nonequilibrium mass transfer process among three phases in soil remediation processes, such as soil vapor extraction (Armstrong et al., 1994; Hoeg et al., 2004) and air sparging (Rahbeh and Mohtar, 2007; VanAntwerp et al., 2008). However, in vapor intrusion models it is assumed that equilibrium exists among three phases. One recent study

considered rate-limited mass transfer between gas and liquid phases for modeling vapor intrusion under the condition of aerobic biodegradation (Chen et al., 2016). In addition, very few physical experiments have been carried out to observe the non-equilibrium phenomena in the column study or the field work, because the precise sampling of contaminants in three phases remains problematic (Cho et al., 1993). This study aimed to formulate a one-dimensional non-equilibrium numerical model with groundwater table fluctuations. The concentrations of VOC contaminants before and after fluctuation commenced were calculated and compared to show the influence of groundwater table fluctuations on the release of VOC contaminants from groundwater to the atmosphere. The difference between the equilibrium and nonequilibrium models was also calculated and compared for different soil characteristics. Finally, the effects of groundwater table fluctuation amplitude, VOC type, soil type, and rainfall events, on the release of contaminants from groundwater to the vadose zone and finally to the crawl space of a building were taken into consideration.

2.The transport model 2.1 The conceptual model In the process of the groundwater table fluctuations, the vertical movement of pore gas and water is dominant, therefore, a one-dimensional (vertical) direction was considered suitable for describing the physical phenomenon. Such model was validated by a soil column experiment, discussed in Section 4. The schematic diagram of the

numerical model is shown in Fig. 1. The position of z = 0 was the top of the saturated zone, while z = L was the ground level. A ventilation opening of 1 m was located above ground to simulate the crawl space above the contaminated groundwater, assuming continuous contaminant ventilation of the opening. The atmospheric pressure fluctuation in the crawl space was ignored in this paper. The thickness of the unsaturated zone including the capillary fringe was 5 m, and the contaminant concentration at the lower boundary was deemed constant to represent the DNAPL contaminated site with semi-infinite groundwater contamination source. In the simulation process, the ventilation opening above the soil surface was simplified as a well-mixed container, and the non-equilibrium mass transfer of contaminants only occurred in the capillary fringe and vadose zone. Most vapor intrusion models consider that contaminants in the gaseous phase, liquid phase and solid phase reach equilibrium in the subsurface environment (Costanza-Robinson et al., 2013; DeVaull, 2007; Luo et al., 2015; Shen et al., 2012b). However, we hypothesize that the mass transfer process may not reach equilibrium among these three phases, and the interfacial mass transfer was considered in this paper. The system was assumed isothermal and heat transfer was not considered. While certain organic contaminants will biodegrade naturally under aerobic conditions in the vadose zone (Davis et al., 2009), chlorinated VOCs typically do not (Kurt and Spain, 2013; USEPA, 2012), and, therefore, biodegradation was not considered in this study.

Fig. 1. The schematic diagram of the numerical model

2.2 Governing water and gas movement equations When the groundwater table fluctuates, water would move in the pore space which could be described by the Richards equation (Richards, 1931):

[Cw (h)  S sw

 w hw  w gks krw hw ]  [ ( +1)]  0  t z w z

where: Cw(h): soil water capacity function, m-1; Ssw: water specific storativity, m-1; θw: volumetric soil moisture content, dimensionless; ϕ: total soil porosity, dimensionless; hw: water pressure head, m;

(1)

t: time, s; z: vertical axis, z = 0 at the top of the saturated zone; ρw: water density, kg/m3; g: the acceleration due to gravity acting in the vertical direction z, m/s2; ks: the intrinsic saturated permeability, m2; krw: the relative water permeability compared to the water saturated permeability, dimensionless; μw: the fluid viscosity, kg/(m·s). When the groundwater table dropped, the bottom of the unsaturated zone would transfer to the saturated zone after imbibition. Therefore, Eq. (1) included the water storativity term, which ensured that the Richards equation remained valid in the saturated zone where Cw(h) was zero (Binning, 1994). The soil moisture content θw and soil water capacity Cw(h) were both determined by soil characteristics, which could be described by the van Genuchten function as follows (van Genuchten, 1980):

w 

   wr   wr [1  ( hw ) n ]m

 mn(   wr ) n (hw ) n 1 , hw  0  Cw (h)   [1  ( hw )n ]m1 0, h  0  w where: α, m, n: van Genuchten parameters, m=1-1/n, dimensionless; θwr: residual water content, dimensionless.

(2)

(3)

The relationship between the relative water permeability and normalized water content was shown as follows (Binning, 1994; Shen et al., 2012b): 1/ m m 2 krw  Se1/2 w [1  (1  Sew ) ]

Sew 

 w   wr    wr

(4) (5)

where Sew was the effective water saturation, dimensionless. Combining Eqs. (1-5), the movement of water could be simulated in the process of groundwater table fluctuations, and the flux of water advection could be calculated by the following equation:

qw   K (h)(

hw  gk k h +1)   w s rw ( w +1) z w z

(6)

where K(h) was the hydraulic conductivity in unsaturated or saturated zones, m/s. Since the soil had relatively high permeability, the variation of gaseous pressure in the unsaturated zone could be neglected during the fluctuations of groundwater table (Guo and Jiao, 2008). Therefore, the pore gas was assumed incompressible, and the continuity equation could be drawn to simulate the gas flow: qw qg  0 z z

(7)

where qg (w) was the gas (water) advection flux, m/s.

2.3 Governing contaminant transport equations in a porous medium For the vadose zone and capillary fringe, the mass balance equations for the transport of contaminants could be described as follows: For the gas phase:

[(   w )Cg ] t



Cg  ( qg C g )  ( Dg )   gw (Cg  Cw H ) z z z

(8)

For the liquid phase:

( wCw )  C (q C ) X  ( Dw w )  w w   gw (Cg  Cw H )   ws (Cw  s ) t z z z Kd

(9)

For the solid phase: [(1   )  s X s ] X   ws (Cw  s ) t Kd

(10)

where: Dg (w): the effective dispersion coefficient of the contaminant in the gas (water) phase, m2/s; Cg (w): the concentration of the contaminant in the gas (water) phase, μg/L; Xs: the concentration of contaminant adsorbed in the soil surface, mg/kg; H: Henry’s law constant, dimensionless; γgw: the mass transfer coefficient between the gas phase and liquid phase, s-1; γws: the mass transfer coefficient between the liquid phase and solid phase, s-1; Kd: the water-solid distribution coefficient, m3/kg; ρs: the soil particle density, kg/m3. In Eqs. (8-9), the effective gas and liquid dispersion coefficients should be determined, which could be calculated by the following equations (Bear, 1972; Millington and Quirk, 1961): 10

Dg  Dg 0

(   w ) 3

2

  | qg |

(11)

Dw  Dw0

10 3 w 2

   | qw | 

(12)

where: Dg0 (w0): the molecular diffusion coefficient of the contaminant in the gas (water) phase, m2/s; ξ: the longitudinal dispersivity or simply dispersivity, m. Our model considered mechanical dispersion because groundwater table fluctuations would induce advective flow for both water and vapors. It should be noted that previous vapor intrusion models at steady state only considered diffusion.

3.Numerical solution 3.1 Model boundary conditions

3.1.1 The boundary conditions for water and gas movement

For the water movement model, when rainfall events were not considered, the top boundary condition could be simplified as a no flux boundary. Therefore, the following equation could be achieved: qw ( z  L, t )  

w gks krw hw h (  1)  0  w |x  L,t 1  0 w z z

(13)

Since the groundwater table would fluctuate at a certain frequency, the bottom boundary could be simplified as the given head boundary, which could be expressed as follows: hw ( z  0, t )  h0

(14)

where h0 was the given hydraulic head, which was also the height of the groundwater table, m. When rainfall at the soil surface was simulated, the following boundary condition could be drawn to consider the effect of rainfall:

qw ( z  L, t )  qrain

(15)

where qrain was the rainfall intensity, m/s. The rainfall event would also cause the groundwater table to rise, but this would be counteracted by the flow of groundwater. Therefore, this effect was not considered in the model. To compare the contaminant distribution with and without groundwater table fluctuations, the groundwater table was kept stable until the distribution of contaminant concentration in three phases reached stability. After that, it started to fluctuate for 5000 h. In order to prevent gas moving across the lower boundary, the lower boundary was set at the lowest position during groundwater table fluctuations. Therefore, the height of the groundwater table h0 could be described by the following equation:

 A, t  0h  h0   2 t  A  A sin( T ), 0h  t  5000h

(16)

where: A: the amplitude of groundwater table fluctuation, m; T: the period of groundwater table fluctuation, h; t: time, h. For the gas flow in the soil, since the saturated zone was located below the lower

boundary during the process of groundwater table fluctuations, there was no gas flow through the lower boundary as described by Eq. (17):

qg ( z  0, t )  0

(17)

3.1.2 The boundary condition for the contaminant transport model

Boundary condition of constant concentration was applied for the lower boundary of the liquid phase as described by Eq. (18):

Cw ( z  0, t )  Cw0

(18)

where Cw0 is the constant concentration in the lower boundary of the model. This lower boundary condition could represent some contaminated sites with DNAPL, which would continue to release contaminants and keep the concentration of groundwater relatively stable. Since there was no gas at the lower boundary of the model, there was no need to define the lower boundary for the contaminant concentration in the gas phase. For the top boundary condition of the contaminant transport model, the gas volume above the soil surface in the crawl space of a building must be considered. We assumed that the contaminant in the gas phase was well mixed in the crawl space, and ventilation was simulated as an airflow blowing in without contaminants. Therefore, the following equation could be drawn to calculate the variation of contaminant concentration in the crawl space: J top Aarea  Qwind Cgout  Vtop

dCgout dt

(19)

where: Jtop: the flux of contaminant through the soil surface to the crawl space, mg/(m2·s); Aarea: the cross-sectional area of the research field, m2; Qwind: the wind blowing flux across the crawl space, m3/s; Cgout: the average contaminant concentration in the crawl space, μg/L; Vtop: the gas volume of crawl space, m3. Jtop was mainly caused by the gas diffusion and advection, which could be expressed by the following equation: J top   Dgtop

Cg z

|z  L  qgtop Cgtop

(20)

where: Dgtop: the gas dispersion coefficient at the soil surface, m2/s; qgtop: the gas flux from the unsaturated zone to the crawl space, m/s. Based on Eq. (10), the adsorbed concentration at the soil surface was totally determined by the aqueous concentration. Therefore, there was no need to give the boundary conditions for the concentration of the contaminant in the solid phase.

3.2 Model initial conditions The initial conditions for the simulation had two parts, the water distribution profiles and the contaminant distribution profiles. At the start of simulation (t = 0 h), the condition for soil water was assumed to be hydrostatic (Shen et al., 2012b), assuming no gas and liquid flow occurred before groundwater table fluctuation. On the other hand, the concentration distribution of contaminant in three phases reached stable

at the start of simulation, and the released flux from unsaturated zone to the crawl space was constant.

3.3 Model discretization In this study, the finite difference method was applied to solve these equations. In the model, the total length of the unsaturated zone was 5 m, which was discretized into 500 grids with the grid length of 1 cm. The lengths of crawl space above soil surface was 1 m (Fig. 1). The time step was 1 h during groundwater table fluctuations in order to catch its movement with the fluctuation period of 12 h. When solving Eqs. (8-9), one difficult problem was to discretize the advection term, and how to choose the concentration in the upstream and downstream sides properly. Many mathematical methods have been developed to solve this problem, such as the upwind scheme, leap-frog scheme and Lax-Friedrichs scheme (Morton and Roe, 2001; Roe, 1989). In this paper, the upwind scheme was applied to determine the concentration of the contaminant in the advection term as shown in Eqs. (21-22):

C t 1

g ( w)i 

C t 1

g ( w)i 

1 2

1 2

Cgt (1w)i , q t 1 1  0 g ( w)i   2   t 1 t 1 Cg ( w)i 1 , qg ( w)i  1  0  2

(21)

Cgt (1w)i 1 , q t 1 1  0 g ( w)i     t 1 t 1 2 Cg ( w)i , qg ( w)i  1  0  2

(22)

where i was the index number for the grid, and (i+1/2) and (i-1/2) meant the position of the grid surface. Eqs. (21-22) illustrates that the concentration at the grid surface was totally determined by the concentration at the upstream grid.

3.4 Mass transfer coefficients among three phases In the contaminant transport model considering non-equilibrium mass transfer, the gas-water and water-solid mass transfer coefficients were two key parameters that determined the equilibrium state among these three phases. Therefore, it was necessary to determine these two parameters based on previous studies. When calculating the gas-water mass transfer coefficient γgw (γgw = Kgwagw), the interfacial sorption coefficient (Kgw) and the gas-water interfacial area (agw) were two key parameters that needed to be determined. Previous studies found that the superficial velocity had significant influence on the interfacial sorption coefficient Kgw, and many empirical equations had been developed based on experimental results (Cho et al., 1993; Gidda et al., 2011). However, in the natural environment, gas and water phases may or may not be static, which was not considered in most equations (Roberts et al., 1985; Shulman et al., 1955). Therefore, the following equation was chosen to calculate Kgw, because it considered both static and non-static conditions (Szatkowski et al., 1995):

Sh a =0.023Sca0.5 +0.849Re0.861 Sc0.5 a a

(23)

where: Sha: Liquid side Sherwood number, dimensionless; Sca: Liquid side Schmidt number, dimensionless; Rea: Liquid side Reynolds number, dimensionless. Eq. (23) ignored the gas side mass transfer resistance because the gas diffusion coefficient was about four orders of magnitude greater than the liquid diffusion

coefficient, inducing a much higher mass transfer coefficient in the gas side (Sherwood and Holloway, 1940; Szatkowski et al., 1995). Therefore, the calculated liquid side gaswater sorption coefficient by Eq. (23) was used as the total gas-water sorption coefficient Kgw. The gas-water interfacial area was another parameter that controlled γgw. Peng and Brusseau (2005) used gas-phase tracer tests to measure the gas-water interfacial area and Eqs. (24-26) were drawn to calculate the gas-water interfacial area, which were applied in this study.

 ( w   m ) 21M  M 6(1   ) agw  {1  [ ] } d50 

(24)

  14.3ln(U )  3.72

(25)

0.098U  1.53,U  3.5 M  1.2,U  3.5

(26)

where: d50: the median soil particle diameter, m; θm: the water content when monolayer water covered solid grains, dimensionless; U: the soil uniformity coefficient, U = d60/d10, dimensionless; β, M: parameters that needed to be determined, dimensionless. The water-solid mass transfer coefficient was correlated with the water-solid distribution coefficient as shown below (Brusseau and Rao, 1991):

 ws  3.89 105  Kd0.53  (1   ) s where γws is in the unit of s-1 and Kd is in the unit of m3/kg .

(27)

3.5 The distribution coefficients among the three phases Previous research found that the water-solid distribution coefficient has a linear relationship with the percentage of soil organic matter (Kile et al., 1995), which could be described by the following equation:

Kd  Koc foc

(28)

where: Koc: the organic carbon-normalized partition coefficient, m3/kg; foc: the organic carbon content of soil, dimensionless. Henry’s law controls the distribution of contaminants between gas and water phases, which is mainly determined by temperature. In this paper the temperature was assumed isothermal at 20 oC, therefore Henry’s law constant (H) was 0.972 for CCl4 according to Chen et al. (2012). It should be regarded that higher temperatures would lead to a higher Henry’s law constant, which would lead to a larger released flux from groundwater to the vadose zone to the crawl space of a building.

3.6 Input parameter values

3.6.1 General simulations

The general base-case scenario input parameter values used are displayed in Table 1. The parameters α, m, n for the van Genuchten equation indicated that the soil was mainly formed by medium sand (Leij et al., 1997). The default fluctuation amplitude of groundwater table was 0.05 m, which could happen in coast aquifers and some inland

area by tidal effect (Marechal et al., 2002; Singaraja et al., 2018). Many parameters in the model are sensitive to temperature, especially the Henry’s law constant (Chen et al., 2012). The parameters in Table. 1 were used in all simulations, unless stated otherwise.

Table 1 Base case input parameter values for general simulations Parameter

Value

Total length of the unsaturated zone and capillary fringe

5m

Total length of the ventilation opening above the soil surface

1m

Cross-sectional area of research field (Aarea)

49 m2

Wind blowing flux in the crawl space (Qwind)

0.0136 m3/s

Temperature

20 oC

Henry’s law constant of CCl4 (H)

0.972

Soil porosity (ϕ)

0.374

Specific storage of groundwater in the saturated zone (Ssw)

10-4 m-1

Residual water content (θwr)

0.055

van Genuchten parameter (n)

3.24

van Genuchten parameter (m)

0.69

van Genuchten parameter (α)

3.32 m-1

Water density (ρw)

103 kg/m3

Acceleration of gravity (g)

9.8 m/s2

Intrinsic saturated permeability (ks)

8.42×10-12 m2

Water dynamic viscosity (μw)

10-3 kg/(m∙s)

Longitudinal dispersivity (ξ)

0.005 m

Gas-water mass transfer coefficient (γgw)

Calculated by Eqs. (25-27)

Median soil particle diameter (d50)

4.0×10-4 m

Water content when monolayer water covered solid grains (θm) 0.00137 Soil uniformity coefficient (U)

5

Organic carbon content of soil (foc)

0.01

Organic carbon-normalized partition coefficient (Koc)

0.0631 m3/kg

Gas molecular diffusion coefficient of CCl4 (Dg0)

8.1×10-6 m2/s

Water molecular diffusion coefficient of CCl4 (Dw0)

8.6×10-10 m2/s

Soil particle density (ρs)

2.55×103 kg/m3

Initial aqueous concentration of CCl4 in groundwater (Cw0)

5 mg/L

Default groundwater fluctuation period (T)

12 h

Default groundwater fluctuation amplitude (A)

0.05 m

3.6.2 VOC type

Four common VOC contaminants of concern, including CCl4, PCE, TCE and benzene, which were selected as target contaminants. The basic parameters of these contaminants are listed in Table A1. When comparing the different types of VOC, the organic content of soil (foc) was all assumed to be 1%. The water-solid distribution

coefficients of different contaminants were calculated by Eq. (28), and the mass transfer coefficients between liquid and solid were calculated by Eq. (27).

3.6.3 Soil types and permeability

Different types of soil may affect the water distribution in soil, which is related to the van Genuchten parameters. Four types of soil were considered, which included coarse sand and gravel, medium sand, silt and clay. The representative parameters for these soils were shown in Table A2.

3.6.4 Effects of infiltration

When considering infiltration, we applied rainfall amounts of 200 mm/y and 500 mm/y, with the rainfall intensity and rainfall times shown in Table A3.

4.Model validation by laboratory experiments In order to check the validity of the model, the numerical model was used to interpret the measured released flux from a previous soil column experiment carried out by Guo (2015). In that experiment, Quikrete® play sand sieved to > 50 mesh size (> 0.3 mm) was placed in a tank with dimensions of 180 cm × 60 cm × 10 cm. The characteristic parameters of the soil, the contaminant (TCE), and column are listed in Table A4. The flow rate of groundwater was 0.32 L/h and the level of groundwater kept increasing until it reached 0.90 m. A period of time was delayed before the groundwater table fluctuation. After that, the groundwater started fluctuation with the help of water

level table elevation control. During the process of groundwater table fluctuations, a 5 L air bag was connected with the contaminant water reservoir to keep the pressure neutral, and a relatively constant concentration of dissolved TCE in the reservoir was kept during the process of groundwater table fluctuations. The TCE concentration in the sweep gas was measured, and the released flux of TCE was calculated during the fluctuations of groundwater table. In the numerical model, the lower boundary for the liquid contaminant was the constant concentration (Cw = 0.9 mg/L), and the upper boundary for the gaseous concentration was also the constant concentration (Cg = 0 mg/L). Model validation was carried out with a fluctuation period of about 310 h and a fluctuation amplitude of about 30 cm. It could be observed that the simulated flux had similar variation trends with the measured flux (Fig. 2a). The flux increased with the decrease of groundwater table, while it decreased with the increase of groundwater table. When the groundwater table dropped, the top part of the saturated zone would transfer to the capillary fringe and unsaturated zone. Because this part had a relatively high TCE concentration, the gaseous concentration of TCE in this part would become very high. The flux would increase due to gas dispersion during the drainage process. However, a small delay existed between the measured and simulated flux, and the simulated flux would drop more slowly than the measured flux. This may be caused by an underestimation of the gaseous dispersion coefficient. In the tank experiment, the gas and water movement was not ideal one-dimensional, and there may be gas and water horizontal movement

both in the saturated zone and capillary fringe (Silliman et al., 2002), which could lead to a transverse dispersivity. To consider the effect of horizontal gas movement, the following equation was drawn to modify gas dispersion coefficient:

 g10/3 Dg  aDg 0 2   | qg | 

(29)

where a was the enhancement factor that represented the effect of horizontal flow on the gas dispersion coefficient. In this experiment, a was chosen as 2.5, and the model validation result is shown in Fig. 2b. It could be observed that the revised model could better simulate the variation of gaseous flux from the soil to the atmosphere. In the real environment, the gas dispersion coefficient may also be increased due to the horizontal gas movement, which should be taken into consideration in future research.

Fig. 2. Comparison of simulated and measured normalized released fluxes of TCE during groundwater table fluctuation experiments conducted by Guo (2015) (a: simulated flux with no revision of the gaseous dispersion coefficient; b: simulated flux with a revised gaseous dispersion coefficient)

5.Results and discussion 5.1 Contaminant concentration distribution in the vadose zone Before the fluctuation of groundwater table commenced, the soil moisture content decreased significantly from the capillary fringe to the top of the unsaturated zone (Fig. 3b). The length of capillary fringe Hc could be calculated by the following equation (Dexter and Bird, 2001): Hc 

1 1 1 m ( )  m

(30)

According to Eq. (30), the length of the capillary fringe was calculated to be 0.34 m, which was relatively short due to the texture of medium sand. On the other hand, before the fluctuation of the groundwater table, the concentration at the top of the capillary fringe was only 0.65% of the equilibrated gaseous CCl4 concentration at the groundwater table (Fig. 3a). The sharp reduction of contaminant concentration through the capillary fringe indicated that the capillary fringe could “block” the transport of volatile organic contaminants from the groundwater table to the unsaturated zone, which was shown by some other studies (Mccarthy and Johnson, 1993; Shen et al., 2013b). The “blocking” effect of the capillary fringe was caused by the variation of moisture content from the capillary fringe to the unsaturated zone (Fig. 3b), which resulted in a significant decrease of the total effective diffusion coefficient and correspondingly a high concentration gradient according to the conservation of mass. However, the “blocking” effect was less significant in silt or clay, because these aquifer materials have a less significant change of moisture content with distance through the

capillary fringe (Shen et al., 2013b). After the groundwater table started fluctuating, the gaseous concentration of CCl4 in the unsaturated zone increased significantly, as the “blocking” effect by the capillary fringe became less significant (Fig. 3a). The concentration at the top of capillary fringe (0.34 m above the initial groundwater table) increased from 0.029 mg/L to 0.29 mg/L at 600 h, and finally reached 0.59 mg/L at 5000 h after fluctuations, thus increasing the contaminant concentration in the unsaturated zone. The flux increase during groundwater table fluctuations was mainly caused by a relatively significant variation in soil moisture content in the capillary fringe (Fig. 3b), which caused an enhanced dispersive flux (Section 5.2). Since the gas-water interfacial area (agw) was strongly correlated with the soil moisture content, it could be observed that the gas-water mass transfer coefficient (γgw) decreased from the top to the bottom of the unsaturated zone (Fig. 3c). The range of γgw was from 1.8×10-4 s-1 to 2.3×10-3 s-1, which was in the same order of magnitude with the water-solid mass transfer coefficient (γws = 1.3×10-3 s-1) in this study. This simulation result indicated that the fluctuations of groundwater table would weaken the “blocking” effect for the reduction of contaminant concentration from groundwater to the unsaturated zone.

Fig. 3. The distribution profiles of gaseous concentration, soil moisture content and gas-water mass transfer coefficient (γgw) before and after groundwater table fluctuations (a: the distribution profiles of gaseous CCl4 concentration in the unsaturated zone; b: the distribution profiles of soil moisture content; c: the distribution profiles of γgw)

The released flux of CCl4 from the unsaturated zone to the atmosphere was 5.3×10-9 g/(m2·s) before groundwater table fluctuation. However, it gradually increased to about 9.3×10-8 g/(m2·s) after fluctuations for 5000 h, which is more than one order of magnitude greater (Fig. 4a). The significant increase in flux to the atmosphere

indicates that groundwater table fluctuations in sandy soil aquifers would greatly enhance the release of VOCs from groundwater to the atmosphere, thus greatly enhancing the potential health risk for humans. The enhancement of VOC emission rate by groundwater table fluctuation was previously observed by Guo et al. (2019), who performed field experiments and found that the TCE emission rate increased from 4 μg/(m2·d) to 8 μg/(m2·d) when the groundwater table fluctuation amplitude increased from 1.5 cm to 10 cm. Illangasekare et al. (2014) also performed soil tank experiments and found that the TCE airflow concentration increased in both water drainage and imbibition processes, which also indicates the enhancement effect of groundwater table fluctuations on VOC flux from groundwater. The flux of CCl4 to the atmosphere did not change significantly at the start of groundwater table fluctuations. It increased gradually at 500 h after fluctuations started (Fig. 4a), because the higher concentration caused by groundwater table fluctuations needed time to transport from the bottom of the unsaturated zone to the top. As shown in Fig. 3a, the gaseous concentrations had increased significantly at the bottom of the unsaturated zone at 600 h, while it had not change significantly at the top of the unsaturated zone. The gaseous concentration at the top of unsaturated zone gradually increased when the fluctuations kept for 2400 h, which shows that much time is needed to increase the gaseous concentration at the top of the unsaturated zone as well as the released flux. On the other hand, the fluctuation amplitude of released flux increased with the continuous fluctuations. After groundwater table fluctuations for 5000 h, the

amplitude of released flux was 8.4×10-9 g/(m2·s), which was 9.1% of the average flux (Fig. 4b). A phase difference of π/6 could be observed between the variation of released flux and the variation of groundwater table. The phase difference indicates that the peak of contaminant concentration in the soil and the released flux to the atmosphere happens when the groundwater table was approaching the highest position. It also indicates that the variation of released VOCs to the atmosphere is both determined by the groundwater table position and its variation rate. However, some experimental results have indicated that the released flux out from the unsaturated zone increases dramatically when the groundwater table suddenly dropped (Guo, 2015; Illangasekare et al., 2014; Mccarthy and Johnson, 1993; Werner and Hohener, 2002). This was mainly caused by the decreased water content in the capillary fringe and top of saturated zone, which increased the exposure of gas phase with water phase with high contaminant concentration. The contaminant in gas and water phases quickly reached equilibrium, resulting in elevated soil gas concentration. Some enclosed contaminant in air bubbles released to the atmosphere during the drainage process may also increase the released flux.

Fig. 4. The variation of released CCl4 flux to the atmosphere before and after groundwater fluctuations (a: time scale -500~5000 h; b: time scale 4950~5000 h)

5.2 The variation of gaseous concentration caused by different mechanisms The variation of gaseous concentration caused by different mechanisms at 2.5 m below the soil surface were calculated and compared to determine the main factors that affect VOC flux during groundwater table fluctuations (Fig. 5). This depth was selected because it was in the middle of unsaturated zone, which is representative for the variation of contaminant concentration at different depths. It could be observed that the variation of gaseous concentration caused by dispersion gradually increased from 3.2×10-9 mg/(L∙s) to 5.4×10-9 mg/(L∙s) when time increased from 250 h to 300 h, which was the main mechanism for the increase of contaminant concentration from groundwater to the atmosphere. The increase of gaseous concentration by dispersion was mainly caused by two factors: (1) The fluctuations of groundwater table increased the flowing flux of gas and water in the unsaturated zone, which would significantly increase the mechanical dispersion of contaminant in both gas and water; (2) When

groundwater table dropped, gas entered the position which was originally water, and the mass transfer from water to gas created a large gaseous concentration in this position. The large contaminant concentration in the bottom of unsaturated zone caused a larger concentration gradient, which finally caused an enhanced dispersive flux. Because the unsaturated zone above the capillary fringe had a larger effective dispersion coefficient compared with the capillary fringe, a larger gaseous concentration in the unsaturated zone would facilitate dispersion due to Fick’s law. This mechanism is revealed by Fig. A1, which shows that the CCl4 flux to the soil surface would also increase significantly after groundwater table fluctuations without considering mechanical dispersion.

Fig. 5. The variation of gaseous concentration caused by different mechanisms during groundwater table fluctuations at 2.5 m below the soil surface (time scale 250~300 h)

The mechanical dispersion of contaminant was strongly correlated with the longitudinal dispersivity (ξ). Sensitivity analysis of the longitudinal dispersivity showed that larger ξ would create higher VOCs flux to the atmosphere, while ignoring

the mechanical dispersion would underestimate the VOCs flux released from the groundwater (Fig. A1). On the other hand, the comparison of the dispersion coefficient caused by diffusion and mechanical dispersion shows that the dispersion in the gas phase was mainly controlled by diffusion (Fig. A2a). The dispersion in the liquid phase was controlled by diffusion in most parts of unsaturated zone, while near the capillary fringe it was mainly controlled by mechanical dispersion (Fig. A2b). This comparison indicates that the mechanical dispersion in the liquid phase near the capillary fringe caused the difference in different dispersivity that was shown in Fig. A1. The mass transfer of contaminants among the three phases plays a negative role on the increase of contaminant flux to the atmosphere (Fig. 5), which was described as a “retardation factor” in previous studies (Brusseau et al., 1997; Costanza-Robinson et al., 2013; Kim et al., 2005). On the other hand, the advective flux would only cause sinusoidal fluctuations of flux in the soil and above the soil surface. Though the transient advective flux could be greater than the diffusive flux, the net contribution of the advective flux was much smaller than the diffusive flux. You and Zhan (2013) developed an equilibrium model and found a similar result. They also concluded that the advective flux contributed comparably with the diffusive flux only when the gasfilled porosity was less than 0.05.

5.3 Comparison between equilibrium and non-equilibrium model The degree of equilibrium between gas and liquid phases (DG-L) was defined by the following equation:

DG  L 

Cw H Cg

(31)

The degree of equilibrium between water and solid phases (DL-S) was defined by the following equation: DL  S 

Xs K d Cw

(32)

The calculated degree of equilibrium between gas and liquid (DG-L) was very close to 1 at all depths in the unsaturated zone, which indicated that the contaminant almost reached equilibrium in most parts of the unsaturated zone (Fig. 6a). Among three depths of 0.005 m, 2.500 m and 4.875 m, DG-L varied most significantly at depth of 0.005 m. When the groundwater table moved upwards across the center position, DG-L reached its minimum of 0.995 at depth of 0.005 m, indicating that the contaminant was transferring from gas to liquid. When the groundwater table moved downwards across the center position, DG-L reached its maximum of 1.005, which indicated that VOCs was transferring from liquid to gas. Similar with DG-L, DL-S was also very close to 1 at all depth of unsaturated zone, which indicated that the contaminant between liquid and solid phases was almost in equilibrium in all parts of the unsaturated zone. The maximal variation of DL-S happened near the soil surface, which was between 0.991 and 1.008 at depth of 0.005 m (Fig. 6b). When the groundwater table rose, the pore water with high contaminant concentration flowed upwards, which caused a higher concentration in the water phase. This process would induce the VOCs to transfer from liquid to solid. Therefore, DL-S decreased with

the rise of the groundwater table and vice versa, as shown in Fig. 6b. From Fig. 6a and Fig. 6b, it can be concluded that slight non-equilibrium mass transfer would occur between gas and water phases, and between water and solid phases. However, this non-equilibrium phenomenon mainly occurs at a short distance above the groundwater table as well as the position below the soil surface. Beyond these regions, the contaminant was almost in equilibrium among three phases. However, in this paper the gas-water mass transfer coefficient and water-solid mass transfer coefficient were both determined by previous experiments (Brusseau and Rao, 1991; Peng and Brusseau, 2005; Szatkowski et al., 1995), which will vary depending on different contaminants and different soil textures. As these two mass transfer coefficients were the key parameters controlling the non-equilibrium mass transfer of VOCs in the unsaturated zone, more research is needed to carefully quantify these two mass transfer coefficients for different contaminants and different soil characteristics.

Fig. 6. The degree of equilibrium among three phases at different depths (a: between gas and liquid phases, b: between liquid and solid phases)

For soils composed with medium sand, the difference in VOC flux between equilibrium and non-equilibrium models was relatively small (Fig. A3). For example, for fine sand with organic content of 1%, the difference of released flux when the fluctuations kept for 5000 h by equilibrium model was about 4.5% lower than that of the non-equilibrium model. Sensitivity analysis of gas-water mass transfer coefficient (γgw) showed that the released flux would not change significantly if γgw decreased to 1% of the original value, while the released flux had a decrease of 33.3% when γgw decreased to 0.1% of the original value (Fig. A4). This result shows that the mass transfer process between gas and liquid would not influence the released flux until it was three orders of magnitude lower than the calculated value. On the other hand, lower water-solid mass transfer coefficients (γws) lead to quicker increases in flux levels during groundwater table fluctuations (Fig. A5), which shows that γws was more sensitive than γgw in a sandy aquifer. The comparison of equilibrium and nonequilibrium models was also carried out in soils composed with silt and clay, for which it was found that equilibrium was also almost fully reached. However, a non-marginal difference between equilibrium and non-equilibrium models were discovered in the case of coarse sand. Results showed that in coarse sand, the released flux calculated by a non-equilibrium model was lower than that of the equilibrium model (Fig. 7a). At 5000 h, the average released flux of CCl4 to the atmosphere calculated by non-equilibrium model was 2.7×10-7 g/(m2·s), which was 18% lower than that calculated by the equilibrium model. This finding indicates that the

traditional equilibrium model might somewhat overestimate the potential human health risk under this condition. During 4950 h to 5000 h without rainfall, the released flux calculated by equilibrium model could reach 2.8×10-7 g/(m2·s) to 3.5×10-7 g/(m2·s), while the released flux was only 2.1×10-7 g/(m2·s) to 3.3×10-7 g/(m2·s) for the nonequilibrium model (Fig. 7b). Though the average flux of the non-equilibrium model was smaller than that calculated by equilibrium model, the fluctuation amplitude of the non-equilibrium model was larger. The discrepancy in coarse sand was mainly caused by a relatively smaller gas-water interfacial area for coarse sand, which reduces the gaswater mass transfer coefficient in the capillary fringe.

Fig. 7. The effect of selecting equilibrium and non-equilibrium models on the released flux of CCl4 to the atmosphere in coarse sand (a: time scale -500~5000 h; b: time scale 4900~5000 h)

5.4 The effect of different parameters on the released flux

5.4.1 Groundwater table fluctuation amplitude

Previous studies have shown that groundwater fluctuation amplitudes can be as much as 1.5~10 cm, especially in areas near the coast (Guo et al., 2019; Marechal et al., 2002; Singaraja et al., 2018). The numerical model showed that the fluctuation amplitude of the groundwater table had strong effect on the released flux of VOCs to the atmosphere (Fig. 8). When the fluctuation amplitude was 0.05 m, the average released flux increased to 9.28×10-8 g/(m2·s) at 5000 h after fluctuations, while it significantly decreased to 5.35×10-8 g/(m2·s) when the fluctuation amplitude decreased to 0.02 m. The simulation result indicated that larger fluctuation amplitude would induce larger fluxes. Though the fluctuation amplitude of flux became larger, the average released flux remained relatively close for the groundwater table fluctuation amplitude of 0.05 m and 0.10 m (Fig. 8b). Because the soil moisture content would not respond to the variation of groundwater table due to a hysteresis effect, the released flux to the atmosphere does not always increase with the increase of fluctuation amplitude. Therefore, careful estimation should be made to identify the potential risk for groundwater table fluctuations at different fluctuation amplitudes.

Fig. 8. The effect of groundwater fluctuation amplitude on released flux of CCl4 to the atmosphere (a: time scale -500~5000 h; b: time scale 4900~5000 h)

5.4.2 Different VOC contaminants

The effects of groundwater fluctuations on the flux of four different VOC types was examined using CCl4, PCE, TCE and benzene as target contaminants (Section 3.6.2). The gas and liquid molecular diffusion coefficients of all four contaminants are similar, so their influences were limited. However, the Henry’s law constant and the water-solid distribution coefficient varied significantly among different contaminants (Table A1). It was found that the average CCl4 flux as well as its fluctuation amplitude was much larger than that of PCE, TCE and benzene after groundwater table fluctuations commenced (Fig. 9). This was mainly caused by the larger Henry’s law constant of CCl4 as well as its small sorption capacity on soil. Although PCE has the second largest Henry’s law constant, the water-solid distribution coefficient of PCE is much larger than the other three contaminants, which caused a larger reduced flux by adsorption on soil surface. It could be concluded that contaminant with larger Henry’s law constant and smaller sorption capacity on soil would have greater released flux after groundwater table fluctuations. It should be considered that, some contaminants may biodegrade whereas other will not, however, this aspect was not considered in the present study.

Fig. 9. The released flux of different contaminants to the atmosphere before and after groundwater table fluctuations (a: time scale -500~5000 h; b: time scale 4900~5000 h)

5.4.3 Different soil types and permeability

The numerical model showed that groundwater table fluctuations increased the release of CCl4 to the atmosphere most significantly in soils composed with coarse sand and gravel (Fig. 10a), for which the flux increased from 7.08×10-9 g/(m2·s) to 2.70×107

g/(m2·s) after fluctuations for 5000 h. For medium sand, the flux increased from

5.27×10-9 g/(m2·s) to 9.28×10-8 g/(m2·s) after fluctuations for 5000 h. It could be deduced that groundwater table fluctuations had higher impact on the increase of released flux for coarse sand and gravel than medium sand. However, for silt and clay, groundwater table fluctuations did not enhance the release of VOCs from groundwater to the atmosphere, and the released flux of VOCs remained nearly constant before and after groundwater table fluctuations (Fig. 10b). This phenomenon is related with the distribution of soil moisture content. In coarse sand and medium sand, the soil moisture content changed significantly in the capillary fringe, and the fluctuations of

groundwater table would cause a significant variation of soil moisture content at the bottom of unsaturated zone (Fig. 3b). The drop of groundwater table would create a zone with high contaminant concentration in the gas phase, which would gradually cause a higher gaseous concentration in the unsaturated zone by dispersion and advection as fluctuations proceeded. However, in silt and clay the variation of soil moisture content during groundwater table fluctuations was insignificant (Fig. A6).

Fig. 10. The released flux of CCl4 to the atmosphere in soils with different textures (a: linear scale (time scale -500~5000 h); b: log scale (time scale 4900~5000 h))

5.4.4 Effects of infiltration

When rainfall events were introduced during the process of groundwater table fluctuations, the released flux to the atmosphere significantly decreased (Fig. A7a). During the rainfall events, the increase of soil moisture content at the top of unsaturated zone blocks the dispersion of contaminants in the gas phase from the unsaturated zone to the atmosphere, which causes a sudden drop of released flux (Shen et al., 2012b).

However, after the rainfall ceases, the flux quite quickly recovers to the level before rainfall. With increased number of rainfall events, the flux drop would naturally become more frequent (Fig. A7b). The released flux was quite close for rainfall intensity of 200 mm/year and 500 mm/year after groundwater table fluctuations for 5000 h. It can be concluded that the influence of rainfall was quite temporary, and could not affect the released flux after the rain stopped for about 100 h.

5.4.5 The multiplying number of released flux by the regression model

As mentioned in the Sections 5.4.1-5.4.3, many parameters affect VOC flux during groundwater table fluctuation. The key parameters that were determined are the soil permeability (i.e., soil type), groundwater table fluctuation amplitude, Henry’s law constant and the organic carbon-normalized partition coefficient (contaminant characteristics). Therefore, these four parameters were used as known parameters to predict the multiplying number for released flux (N) after 5000 h of groundwater table fluctuation. The following equation was obtained by the multi-parameter regression model to predict the multiplying number (N):

N  124.9477  12.0162lg ks  16.5821H  50.1644Koc  231.8718 A

(33)

The p value of the regression model is 0.0020 (<0.05), which indicates that it effectively predicts N with small error.

6.Conclusions and limitations This study developed a numerical model to study the effect of groundwater table

fluctuations on the transport of VOCs from groundwater to the atmosphere in a sandy aquifer. In this numerical model, the mass transfer processes between gas and liquid phases, and between liquid and solid phases were considered. The model was validated by the groundwater table fluctuation experiment by Guo (2015), and results showed that the numerical model could simulate the variation trend of released flux during groundwater table fluctuations well. However, enhanced gas dispersion was found during the experiment, which might be caused by the horizontal gas movement in the tank. In a sandy aquifer, when VOCs transported from the contaminated groundwater to the atmosphere without the fluctuation of groundwater table, its concentration decreased by more than two order of magnitude from the bottom to the top of the capillary fringe due to the slow diffusion rate in the pore water. However, groundwater table fluctuations eliminate this “blocking” effect and greatly enhance VOC flux to the atmosphere, therefore significantly increasing the potential health risk for human. This enhancement of contaminant emission was mainly caused by enhanced dispersion, while the mass transfer among three phases would weaken this effect. The dispersion in the gas phase was mainly controlled by diffusion, while the mechanical dispersion was dominant near the capillary fringe in the liquid phase. In the sandy aquifer, equilibrium among three phases was almost reached in the unsaturated zone during the groundwater table fluctuations. In coarse sand, a reduced interfacial area between gas phase and liquid phase leads to a relatively small mass

transfer coefficient, which causes a non-marginal overestimation of the released flux by the equilibrium models. Further research is needed to carefully determine the mass transfer coefficient among three phases for different contaminants under different conditions e.g., changes in barometric pressure or the amount of soil organic matter present. The released flux increased significantly when the fluctuation amplitude increased from 0.02 m to 0.05 m. However, a hysteresis effect in soil moisture content limits the increase in VOC flux with larger fluctuation amplitudes. VOCs with higher Henry’s law constants and lower water-solid distribution coefficients will have greater flux levels during groundwater table fluctuations. Groundwater table fluctuations greatly enhance the released flux of VOCs in coarse sand and medium sand. However, they hardly increase the released flux in silt or clay, because the soil moisture content changes little during groundwater table fluctuations. Rainfall events cause a sudden drop in released flux, while its effect is temporal and does not greatly affect the overall increase of released flux during groundwater table fluctuations. Finally, a regression equation was carried out to predict the variation of released flux with groundwater table fluctuation. It should be noted that in the present study the soil (unsaturated zone) was all composed with materials with no soil horizons (layers) or soil heterogeneity, which is unrepresentative of natural soil. Moreover, the adsorption processes (solid-liquid, gasliquid) were assumed to be linear, and the mixing of salt water and fresh water near the

coast, the atmospheric pressure fluctuations, the preferential flow path, the biodegradation effect and some other factors were not considered in this paper. These limitation factors could be investigated in further studies based on this model. The numerical model presented in this paper is the first step toward improved estimation of vapor intrusion risk at contaminated sites with groundwater table fluctuations and mass transfer processes among three phases.

Declaration of Competing Interest The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.

Acknowledgment This work was supported by the National Natural Science Foundation of China (Grant No. 41671316). Moreover, this study was partly funded by China’s National Water Pollution Control and Treatment Science and Technology Major Project (Grant No. 2018ZX07109-003), and the National Key Research and Development Program of China (Grant No. 2018YFC1801300). The data necessary to reproduce the conclusions of

our

manuscript

are

available

on

the

following

site:

https://figshare.com/articles/Influence_of_groundwater_table_fluctuation_on_the_tra nsport_of_volatile_organic_contaminants_in_the_vadose_zone/8856890

References: Armstrong, J.E., Frind, E.O., Mcclellan, R.D., 1994. Nonequilibrium MassTransfer between the Vapor, Aqueous, and Solid-Phases in Unsaturated Soils during Vapor

Extraction.

Water

Resources

Research,

30(2):

355-368.

DOI:Doi

10.1029/93wr02481 Bear, J., 1972. Dynamics of Fluids in Porous Media. Dover Publications, New York. Binning, P.J., 1994. Modeling unsaturated zone flow and contaminant transport in the air and water phases, Princeton University, Princeton, New Jersey. Brusseau, M.L., Popovicova, J., Silva, J.A.K., 1997. Characterizing gas-water interfacial and bulk water partitioning for gas phase transport of organic contaminants in unsaturated porous media. Environmental Science & Technology, 31(6): 1645-1649. DOI:DOI 10.1021/es960475j Brusseau, M.L., Rao, P.S.C., 1991. Influence of Sorbate Structure on Nonequilibrium

Sorption

of

Organic-Compounds.

Environmental

Science

&

Technology, 25(8): 1501-1506. DOI:DOI 10.1021/es00020a022 Chen, F., Freedman, D.L., Falta, R.W., Murdoch, L.C., 2012. Henry's law constants of chlorinated solvents at elevated temperatures. Chemosphere, 86(2): 15665. DOI:10.1016/j.chemosphere.2011.10.004 Chen, Y., Hou, D., Lu, C., Spain, J.C., Luo, J., 2016. Effects of Rate-Limited Mass Transfer on Modeling Vapor Intrusion with Aerobic Biodegradation. Environ Sci

Technol, 50(17): 9400-6. DOI:10.1021/acs.est.6b01840 Cho, H.J., Jaffe, P.R., Smith, J.A., 1993. Simulating the Volatilization of Solvents in Unsaturated Soils during Laboratory and Field Infiltration Experiments. Water Resources Research, 29(10): 3329-3342. DOI:Doi 10.1029/93wr01414 Costanza-Robinson, M.S., Carlson, T.D., Brusseau, M.L., 2013. Vapor-phase transport of trichloroethene in an intermediate-scale vadose-zone system: retention processes

and

tracer-based

prediction.

J

Contam

Hydrol,

145:

82-9.

DOI:10.1016/j.jconhyd.2012.12.004 Davis, G.B., Patterson, B.M., Trefry, M.G., 2009. Evidence for Instantaneous Oxygen-Limited Biodegradation of Petroleum Hydrocarbon Vapors in the Subsurface. Ground Water Monitoring and Remediation, 29(1): 126-137. DOI:10.1111/j.17456592.2008.01221.x DeVaull, G.E., 2007. Indoor vapor intrusion with oxygen-limited biodegradation for a

subsurface

gasoline

source.

Environ

Sci

Technol,

41(9):

3241-8.

DOI:10.1021/es060672a Dexter, A.R., Bird, N.R.A., 2001. Methods for predicting the optimum and the range of soil water contents for tillage based on the water retention curve. Soil & Tillage Research, 57(4): 203-212. DOI:Doi 10.1016/S0167-1987(00)00154-9 Garg, S. et al., 2017. Overview of Natural Source Zone Depletion: Processes, Controlling Factors, and Composition Change. Ground Water Monitoring and Remediation, 37(3): 62-81. DOI:10.1111/gwmr.12219

Gidda, T., Stiver, W.H., Zytner, R.G., 2011. Measuring Aqueous-Air and SorbedAqueous Mass Transfer Coefficients for Application in Soil Vapor Extraction. Journal of Environmental Engineering-Asce, 137(10): 880-888. DOI:10.1061/(Asce)Ee.19437870.0000405 Guo, H.P., Jiao, J.J., 2008. Numerical study of airflow in the unsaturated zone induced by sea tides. Water Resources Research, 44(6): 13. DOI:Artn W06402 10.1029/2007wr006532 Guo, Y.M., 2015. Vapor intrusion at a site with an alternative pathway and a fluctuating groundwater table, Arizona State University, Phoenix, Arizona State. Guo, Y.M., Holton, C., Luo, H., Dahlen, P., Johnson, P.C., 2019. Influence of Fluctuating Groundwater Table on Volatile Organic Chemical Emission Flux at a Dissolved Chlorinated-Solvent Plume Site. Ground Water Monitoring and Remediation, 39(2): 43-52. DOI:10.1111/gwmr.12322 Hoeg, S., Scholer, H.F., Warnatz, J., 2004. Assessment of interfacial mass transfer in water-unsaturated soils during vapor extraction. J Contam Hydrol, 74(1-4): 163-95. DOI:10.1016/j.jconhyd.2004.02.010 Hou, D.Y., Li, F.S., 2017. Complexities Surrounding China's Soil Action Plan. Land Degradation & Development, 28(7): 2315-2320. DOI:10.1002/ldr.2741 Illangasekare, T. et al., 2014. Vapor intrusion from entrapped NAPL sources and groundwater plumes: Process understanding and improved modeling tools for pathway assessment, Colorado School of Mines, Colorado.

Kile, D.E., Chiou, C.T., Zhou, H., Li, H., Xu, O., 1995. Partition of nonpolar organic pollutants from water to soil and sediment organic matters. Environ Sci Technol, 29(5): 1401-6. DOI:10.1021/es00005a037 Kim, H., Lee, S., Moon, J.W., Rao, P.S.C., 2005. Gas transport of volatile organic compounds in unsaturated soils: Quantitative analysis of retardation processes. Soil Science Society of America Journal, 69(4): 990-995. DOI:10.2136/sssaj2003.0208 Kurt, Z., Spain, J.C., 2013. Biodegradation of chlorobenzene, 1,2-dichlorobenzene, and 1,4-dichlorobenzene in the vadose zone. Environ Sci Technol, 47(13): 6846-54. DOI:10.1021/es3049465 Leij, F.J., Russell, W.B., Lesch, S.M., 1997. Closed-form expressions for water retention

and

conductivity

data.

Ground

Water,

35(5):

848-858.

DOI:DOI

10.1111/j.1745-6584.1997.tb00153.x Luo, J., Kurt, Z., Hou, D., Spain, J.C., 2015. Modeling aerobic biodegradation in the capillary fringe. Environ Sci Technol, 49(3): 1501-10. DOI:10.1021/es503086p Marechal, J.C., Sarma, M.P., Ahmed, S., Lachassagne, P., 2002. Establishment of earth tide effect on water-level fluctuations in an unconfined hard rock aquifer using spectral analysis. Curr. Sci., 83(1): 61-64. Mccarthy, K.A., Johnson, R.L., 1993. Transport of Volatile Organic-Compounds across the Capillary-Fringe. Water Resources Research, 29(6): 1675-1683. DOI:Doi 10.1029/93wr00098 Millington, R., Quirk, J.P., 1961. Permeability of Porous Solids. Transactions of the

Faraday Society, 57(8): 1200-&. DOI:DOI 10.1039/tf9615701200 Morton, K.W., Roe, P.L., 2001. Vorticity-preserving Lax-Wendroff-type schemes for the system wave equation. Siam Journal on Scientific Computing, 23(1): 170-192. DOI:Doi 10.1137/S106482759935914x Peng, S., Brusseau, M.L., 2005. Impact of soil texture on air-water interfacial areas in unsaturated sandy porous media. Water Resources Research, 41(3): 8. DOI:Artn W03021 10.1029/2004wr003233 Picone, S., Valstar, J., van Gaans, P., Grotenhuis, T., Rijnaarts, H., 2012. Sensitivity analysis on parameters and processes affecting vapor intrusion risk. Environ Toxicol Chem, 31(5): 1042-52. DOI:10.1002/etc.1798 Rahbeh, M.E., Mohtar, R.H., 2007. Application of multiphase transport models to field remediation by air sparging and soil vapor extraction. J Hazard Mater, 143(1-2): 156-70. DOI:10.1016/j.jhazmat.2006.09.098 Raubenheimer, B., Guza, R.T., Elgar, S., 1999. Tidal water table fluctuations in a sandy ocean beach. Water Resources Research, 35(8): 2313-2320. DOI:Doi 10.1029/1999wr900105 Recknagel, R.O., 1967. Carbon tetrachloride hepatotoxicity. Pharmacol Rev, 19(2): 145-208. Richards, L.A., 1931. Capillary conduction of liquids through porous mediums. Physics-J Gen. Appl. Phys., 1(1): 318-333. DOI:Doi 10.1063/1.1745010

Rivett, M.O., Wealthall, G.P., Dearden, R.A., McAlary, T.A., 2011. Review of unsaturated-zone transport and attenuation of volatile organic compound (VOC) plumes leached from shallow source zones. J Contam Hydrol, 123(3-4): 130-56. DOI:10.1016/j.jconhyd.2010.12.013 Roberts, P.V., Hopkins, G.D., Munz, C., Riojas, A.H., 1985. Evaluating 2Resistance Models for Air Stripping of Volatile Organic Contaminants in a Countercurrent, Packed-Column. Environmental Science & Technology, 19(2): 164173. DOI:DOI 10.1021/es00132a010 Roe, P., 1989. A survey of upwind differencing techniques, 11th International Conference on Numerical Methods in Fluid Dynamics. Springer, pp. 69-78. Shen, R., Pennell, K.G., Suuberg, E.M., 2012a. A numerical investigation of vapor intrusion--the dynamic response of contaminant vapors to rainfall events. Sci Total Environ, 437: 110-20. DOI:10.1016/j.scitotenv.2012.07.054 Shen, R., Pennell, K.G., Suuberg, E.M., 2012b. A numerical investigation of vapor intrusion - The dynamic response of contaminant vapors to rainfall events. Sci. Total Environ., 437: 110-120. DOI:10.1016/j.scitotenv.2012.07.054 Shen, R., Pennell, K.G., Suuberg, E.M., 2013a. Influence of Soil Moisture on Soil Gas Vapor Concentration for Vapor Intrusion. Environ. Eng. Sci., 30(10): 628-637. DOI:10.1089/ees.2013.0133 Shen, R., Pennell, K.G., Suuberg, E.M., 2013b. Influence of Soil Moisture on Soil Gas Vapor Concentration for Vapor Intrusion. Environ Eng Sci, 30(10): 628-637.

DOI:10.1089/ees.2013.0133 Shen, R., Pennell, K.G., Suuberg, E.M., 2014. Analytical modeling of the subsurface volatile organic vapor concentration in vapor intrusion. Chemosphere, 95: 140-9. DOI:10.1016/j.chemosphere.2013.08.051 Shen, R., Suuberg, E.M., 2016. Impacts of Changes of Indoor Air Pressure and Air Exchange Rate in Vapor Intrusion Scenarios. Build Environ, 96: 178-187. DOI:10.1016/j.buildenv.2015.11.015 Shen, R., Yao, Y., Pennell, K.G., Suuberg, E.M., 2013c. Modeling quantification of the influence of soil moisture on subslab vapor concentration. Environ Sci Process Impacts, 15(7): 1444-51. DOI:10.1039/c3em00225j Sherwood, T.K., Holloway, F.A.L., 1940. Performance or packed towers - Liquid film data for several packings. Transactions of the American Institute of Chemical Engineers, 36(1): 0039-0070. Shoushtari, S.M.H.J., Cartwright, N., Nielsen, P., Perrochet, P., 2016. The effects of oscillation period on groundwater wave dispersion in a sandy unconfined aquifer: Sand flume experiments and modelling. Journal of Hydrology, 533: 412-420. DOI:10.1016/j.jhydrol.2015.12.032 Shulman, H.L., Ullrich, C.F., Proulx, A.Z., Zimmerman, J.O., 1955. Performance of Packed Columns .2. Wetted and Effective-Interfacial Areas, Gas-Phase and LiquidPhase Mass Transfer Rates. Aiche J., 1(2): 253-258. DOI:DOI 10.1002/aic.690010220 Silliman, S.E., Berkowitz, B., Simunek, J., van Genuchten, M.T., 2002. Fluid flow

and solute migration within the capillary fringe. Ground Water, 40(1): 76-84. DOI:DOI 10.1111/j.1745-6584.2002.tb02493.x Singaraja, C., Chidambaram, S., Jacob, N., 2018. A study on the influence of tides on the water table conditions of the shallow coastal aquifers. Applied Water Science, 8(1): 11. DOI:UNSP 11 10.1007/s13201-018-0654-5 Szatkowski, A., Imhoff, P.T., Miller, C.T., 1995. Development of a Correlation for Aqueous-Vapor Phase Mass-Transfer in Porous-Media. Journal of Contaminant Hydrology, 18(1): 85-106. DOI:Doi 10.1016/0169-7722(94)00035-G Thomson, N.R., Sykes, J.F., Van Vliet, D., 1997. A numerical investigation into factors affecting gas and aqueous phase plumes in the subsurface. Journal of Contaminant Hydrology, 28(1-2): 39-70. DOI:Doi 10.1016/S0169-7722(96)00044-7 USEPA, 2012. Petroleum hydrocarbon and chlorinated solvents differ in their potential for vapor intrusion, Washington, D.C., pp. 1-13. van Genuchten, M.T., 1980. A Closed-form Equation for Predicting the Hydraulic Conductivity of Unsaturated Soils1. Soil Science Society of America Journal, 44(5): 892-898. DOI:10.2136/sssaj1980.03615995004400050002x VanAntwerp, D.J., Falta, R.W., Gierke, J.S., 2008. Numerical simulation of fieldscale contaminant mass transfer during air sparging. Vadose Zone Journal, 7(1): 294304. DOI:10.2136/vzj2006.0126 Verginelli, I., Baciocchi, R., 2011. Modeling of vapor intrusion from hydrocarbon-

contaminated sources accounting for aerobic and anaerobic biodegradation. J Contam Hydrol, 126(3-4): 167-80. DOI:10.1016/j.jconhyd.2011.08.010 Werner, D., Hohener, P., 2002. The influence of water table fluctuations on the volatilization of contaminants from groundwater. In: Thornton, S.F., Oswald, S.E. (Eds.), Groundwater Quality: Natural and Enhanced Restoration of Groundwater Pollution. Iahs Publication. Int Assoc Hydrological Sciences, Wallingford, pp. 213-218. Yao, Y., Verginelli, I., Suuberg, E.M., Eklund, B., 2018. Examining the use of USEPA's Generic Attenuation Factor in determining groundwater screening levels for vapor

intrusion.

Ground

Water

Monit

Remediat,

38(2):

79-89.

DOI:10.1111/gwmr.12276 You, K.H., Zhan, H.B., 2013. Comparisons of diffusive and advective fluxes of gas phase volatile organic compounds (VOCs) in unsaturated zones under natural conditions.

Advances

in

DOI:10.1016/j.advwatres.2012.11.021

Water

Resources,

52:

221-231.

Highlights 

Groundwater table fluctuation increases the released flux of VOCs from groundwater



Equilibrium models conservatively overestimate the released flux in coarse sand



The increased flux was most significant in sandy aquifer



VOCs of higher volatility and less sorption capacity have larger released fluxes