## Abstract

To quantify the vertical diffusive mass transfer of trichloroethylene (TCE) in a low water-saturated porous medium, diffusion experiments were conducted in soil columns using a medium-size quartz sand at two different water saturations; namely, completely dry and at irreducible water saturation. Advective transport of TCE vapors was mainly generated by density gradients in the soil gas during diffusion of TCE vapors in the vertical direction. The maximum stationary TCE vapor concentrations measured at the column outlet attained only 70 to 75% of the prescribed inlet vapor concentration, indicating a strong influence of density-driven advective mass transport on vertical diffusion of TCE. Without taking into account the density of the vapors, the measured apparent TCE diffusion coefficient was about seven and five times smaller than the coefficients predicted by the model for dry sand and sand at irreducible water saturation, respectively. A numerical transport model was used to interpret the experimental results. The numerical results fit the experimentally measured TCE breakthrough curves well. A sensitivity study demonstrated that the simulated vapor breakthrough curves were not very sensitive to variations in the dimensionless Henry's constant, the gas–solid or water–solid phase partitioning coefficients but were sensitive to modification of the intrinsic permeability and the gas tortuosity–connectivity parameter. Two series of numerical simulations were also conducted to analyze the steady-state conditions of vertical vapor transport. The results highlighted how both the relative density of TCE vapors and the height of the soil column, expressed as a dimensionless Rayleigh number, might reduce the vertical diffusive mass flux.

Vertical mass transfer of trichloroethylene (TCE) in a low water-saturated porous medium was studied in laboratory columns. A numerical transport model was used to interpret the measured vapor breakthrough curves, to determine apparent TCE diffusion coefficients, and to analyze the steady-state conditions of vertical vapor mass transport.

**Industrialization and technological development** during the 19th and 20th century led to the use of dense nonaqueous-phase liquids (DNAPLs) such as chlorinated solvents, e.g., TCE. At locations where these volatile organic compounds (VOCs) were accidentally spilled on the soil during transport or leaked from their storage places, VOCs could migrate vertically (depending on gravity and capillary forces) through the unsaturated zone toward the underlying groundwater. Once a DNAPL has penetrated the aquifer, a certain quantity of pollutant could become trapped as “blobs” or form pools over low-permeability geologic structures. As a result of the high volatility and generally low water solubility of DNAPLs, these blobs and pools can become sources for large DNAPL vapor plumes and long-term solute plumes in the groundwater (Hounslow, 1995; Pankow and Cherry, 1996). Many recent studies have shown the high toxicity of VOCs, such as TCE, and the fact that their particular physical and chemical characteristics can cause a high environmental risk (Brüning and Bolt, 2000; Mensing et al., 2002; Lash et al., 2005). Understanding the mechanisms for DNAPL vapor transport at VOC-contaminated sites is important for reducing such environmental risks but still remains a challenge to date (Shapiro et al., 2004).

In natural soils, advective transport of DNAPL vapors might be generated by density gradients that exist within and along the fringe of the vapor plume (Sleep and Sykes, 1989; Lenhard et al., 1995; Pankow and Cherry, 1996) or by pressure gradients. These pressure gradients may be due, for example, to variations in the atmospheric pressure (Massmann and Farrier, 1992; Pasteris et al., 2002) or to volatilization of the DNAPL (Baehr and Bruell, 1990; Altevogt et al., 2003). Mendoza and Frind (1990) developed a sensitivity analysis of their numerical model to predict the migration of DNAPL in the vadose zone. They showed that density-driven advection is not negligible in the case of an air–compound mixture that has a density 1.15 times higher than that of uncontaminated soil air. Johnson and Kreamer (1994) emphasized that the density limit for notable density-driven advection is about 1.2. This ratio is 1.27 at 20°C for saturated TCE vapor and thus density-driven advective mass transport might be significant for this compound.

Density-driven advection is not always considered in the interpretation of TCE diffusion experiments (Batterman et al., 1995, 1996). Although advection might be a significant transport mechanism for vapors in the unsaturated zone, diffusion remains predominant (Pankow and Cherry, 1996; Choi et al., 2002; Jellali et al., 2003; Bohy et al., 2006; Cotel, 2008). Any of the following types of VOC mass transfer between the solid, water, and gas phases might significantly delay vapor transport in the vadose zone: (i) sorption on porous media, which is characterized by a compound partitioning coefficient; (ii) dissolution of the DNAPL vapor in water, which is characterized by the Henry's law constant; and (iii) retention of the VOC at the soil air–soil water interface, which is characterized by an adsorption coefficient (Conklin et al., 1995; Brusseau et al., 1997; Grathwohl, 1998; Kim et al., 2001, 2005). Kim et al. (2001) showed that all of the grain surfaces for a given sand with a water saturation of 3% were covered with a continuous film of water that prevented any contact between the TCE vapors and the solid matrix, suggesting that, in a porous medium at irreducible water saturation, direct vapor adsorption on a porous medium is infeasible, as further supported by others (Dorris and Gray, 1981; Ong and Lion, 1991; Unger et al., 1996; Goss and Schwarzenbach, 2002). Moreover, retention of VOC at the interface between the soil air and soil water is always weak for faintly polar compounds such as TCE (Hoff et al., 1993) unless the porous medium contains very little water (Kim et al., 2001; Cotel, 2008).

The parameters that influence vapor diffusion are water saturation (Millington and Quirk, 1961; Moldrup et al., 2000), pressure (Altevogt et al., 2003), and temperature (Lyman et al., 1990; Grathwohl, 1998). Among these parameters, water saturation is considered to be the most important because it varies considerably in soils and directly influences the tortuosity of the porous medium and thus the diffusion coefficient of the gaseous phase. Experimental laboratory studies of DNAPL vapor migration have been conducted, but the results of these studies were difficult to interpret because the boundary or initial conditions were not adequately controlled and the different mass transfer mechanisms were not precisely quantified (Batterman et al., 1995, 1996; Altevogt et al., 2003; Emonet, 2003). Furthermore, very few studies have run replicate experiments; doing so would reinforce the quality of the experimental results and improve their interpretation.

This study focused on (i) experiments on the influence of density-driven advection on the vertical diffusive transport of TCE vapors under well-known boundary conditions in a laboratory column using a medium quartz sand at two different water saturations, completely dry and irreducible water saturation, (ii) numerical modeling of the TCE vapor breakthrough curves measured at the column outlet with a coupled flow and transport model for the gas phase to evaluate both the validity of the chosen numerical approach and the influence of the vapor density on the “apparent gas diffusion” coefficient when neglecting the advection term, and (iii) sensitivity studies to analyze, in detail, the influence of selected model parameters on the simulated breakthrough curves and to generalize the experimentally observed vertical steady-state TCE vapor flux by varying both the imposed upstream vapor concentration and the column height.

## Materials and Methods

### Experimental Setup

Controlled diffusion experiments were performed in a laboratory column using medium quartz sand at two different water saturations, completely dry and irreducible water saturation. Several replicates of the experiments were performed and special care was taken to choose adequate initial and boundary conditions of the gas pressure, which may reduce the effect of advective transport of TCE vapors due to pressure gradients, to maintain well-defined boundary conditions of the vapor concentration, to control the parameters that influence the vapor transport, and to precisely measure the transient TCE vapor concentrations.

The diffusion column used in the experiments consisted of a 22-cm-long glass column with an internal diameter of 7 cm filled with a medium-size homogeneous quartz sand (Fig. 1 ). The sand had a low fraction of organic matter (*f*_{oc} = 0.09%, based on the French standard NF X 31–109 [Association Française de Normalisation, 1993]) and a hydraulic conductivity of about 8 × 10^{−4} m s^{−1}, corresponding to an intrinsic permeability of 8.6 × 10^{−11} m^{2} (Table 1 ). This sand was already the subject of several former research studies, and its characteristics were described by Benremita (2002), Jellali et al. (2003), Bohy (2003), and Bohy et al. (2006). The aqueous–solid distribution coefficient and longitudinal dispersivity were quantified in former studies; all other sand properties were measured on the given experimental setup and in preliminary column studies using the same filling procedure (Cotel, 2008).

The porous medium was supported at its lower end by a 0.7-cm-thick sintered glass. Two cavities were created at each end of the column (Fig. 1). In the so-called upstream cavity, the DNAPL was introduced at the beginning of the experiment. The downstream cavity was used to take small gas samples to measure the accumulated vapor mass as a function of time at the outlet section without significantly disturbing the diffusion process within the permeable sand. The parts of the three glass modules of the device were assembled, and their seal was ensured by the use of vacuum grease and two locking rings.

Each glass module included openings to take gas samples, introduce measuring devices, or add the DNAPL. Small glass tubes with internal diameters of 8 mm were equipped with threads and closed by a set of Teflon septa, sealing rings, and caps. These septa were the only materials in contact with the porous medium and fluids.

Unlike the classical diffusion cells that were used, for example, by Alzaydi et al. (1978), Fuentes et al. (1991), Thoma et al. (1992), Batterman et al. (1995), and Solcova and Schneider (2003), the chosen device did not include air circulation in both cavities. Forced air circulation might be useful for homogenizing the vapor concentrations but it introduces local pressure gradients in the outlet chamber and disturbs the progress of the vapor front (Cotel, 2008).

The DNAPL used in the diffusion experiments was TCE. Table 2 summarizes the physical and chemical properties of TCE. This pollutant was chosen because it is the solvent that is found most frequently in groundwater throughout the world (Lerner et al., 1991; Schiedek et al., 1997). Furthermore, TCE is very volatile, fairly soluble, dense, and weakly adsorbed on the surface of mineral grains. Biodegradation of TCE, as reported in the literature, did not play a role in the conducted diffusion experiments because TCE breakdown products, such as *cis*-1,2 dichloroethene, were not detected in the analysis of the gas samples. The physicochemical characteristics of TCE (Table 2) are representative of the DNAPL family.

The TCE was injected as a liquid phase in the upstream cavity. The injected volume was chosen to minimize the air volume between the sintered glass and the DNAPL. Stainless steel needles with an internal diameter of 2 mm were placed in the septa of the two cavities to avoid a pressure rise of the air mixture in the upstream cavity because of the volatilization of the DNAPL (see, for example, Exp. A, Table 3 ). The needles also maintained the gaseous phase at a pressure equal to the ambient air pressure.

^{241}Am, and the detector was NaI. The water saturations were obtained using the modified Beer–Lambert law (Gerard-Marchant, 1998) as follows:

*S*

_{w}(

*z*) (dimensionless) is the water saturation at height

*z*, ε (dimensionless) is the porosity of the porous medium, μ

_{w}′ (cm

^{2}g

^{−1}) is the mean attenuation mass factor of the water, ρ

_{w}(g cm

^{−3}) is the water density,

*d*(cm) is the internal diameter of the glass column, and

*N*

_{s}(

*z*) and

*N*(

*z*) are the number of counts monitored in the dry sand fill and in the wet sand fill at height

*z*.

Note that the fraction of the original photons that is emitted from the γ source, passes through the soil column, and arrives at the detector depends on the distance passed through the glass column and the porous medium and their mass absorption coefficients. The number of photons arriving at the detector is usually called the number of counts. In the present study, a period of 60 s was used to determine the number of counts for each point measurement.

The filling of the column was performed with dried sand added successively in layers of several centimeters and mixed from bottom to top with a rigid thin rod. This method was previously demonstrated by Szenknect (2003) to give packed columns with highly reproducible and uniform physical properties. Gamma densitometry was used to verify that the sand fill of the column was homogeneous.

To obtain uniform water saturation in the sand, the method of Kim et al. (2001) was used. Water (25 g) was injected into a total mass of 670 g of the sand fill. This ratio of water to sand corresponds to an irreducible water saturation of 13%. The water was injected on the top of the sealed column using a syringe. Following the procedure of Kim et al. (2001), the whole apparatus was then placed in an oven and heated at 95°C for about 48 h to allow the water vapor to distribute uniformly along the column. The setup was then slowly cooled to the ambient room temperature of 20°C to avoid water condensation on the glass column. Successive weighing of the whole setup showed that there were no significant losses of water during the heating and cooling phases.

During the diffusion experiments, gas samples were taken both from the downstream cavity (at *z* = 0.13 m) and from the upstream cavity (at *z* = −0.02 m) at increasing time steps. The number of vapor samples for each experiment was chosen to find a compromise between getting fine data resolution of the concentration breakthrough curve at the downstream cavity and creating the least possible disruption of the diffusion process. As an upper limit, a total gas volume sampled during the experiments was fixed at 3% of the air content in the device. For estimated TCE concentrations ranging from 1 to 419 mg L^{−1}, vapor samples of 0.25 mL were taken with a microsyringe. In the case of lower TCE concentrations ranging from 10 μg L^{−1} to 1 mg L^{−1}, sample volumes of 0.5 mL were chosen. It should be emphasized here that the low TCE concentration range is at the higher end of what is generally measured in soil gas samples collected in the field.

Each vapor sample was dissolved in a small glass bottle containing 10 mL of hexane. After sampling, the microsyringe needle was rinsed about 10 times with hexane, and then the bottles were shaken for 1 h. The TCE vapor concentration was determined using a gas chromatography device equipped with an electron capture detector (Clarus 500, PerkinElmer, Waltham, MA; capillary column: Elite 5). The analysis lasted 8.5 min, and a volume of 0.5 μL was injected with a 1/15 split. The temperature of the injector was 150°C, and the oven temperature varied between 40 and 80°C. The temperature of the detector was 300°C. The reproducibility of the analysis had an average standard deviation of 1.9% obtained on the basis of 12 identical gas samples of each standard.

### Mathematical Model Used

To model the TCE vapor transport observed in the experimental setup, a mathematical approach was used that was based on a coupled flow and transport model for the gas phase in the given three-phase (water–gas–porous medium) system.

*z*direction for a partially saturated porous medium at the scale of a representative elementary volume (REV). The porous medium was assumed to be incompressible, and the air content was assumed to be constant with time. Without addition or loss of gas within the REV, the generalized Darcy's law (Bear, 1972) expresses the specific discharge of the gas phase. The flow equation for density-dependent gas flow can be written as

*k*

_{ra}(dimensionless) is the relative gas permeability,

*k**(m

^{2}) is the intrinsic permeability, ρ

_{a}(kg m

^{−3}) is the gas density, μ

_{a}(Pa s) is the dynamic gas viscosity,

*P*

_{a}(Pa) is the gas pressure,

*g*(m s

^{−2}) is the gravity acceleration constant, and

*S*

_{a}(dimensionless) is the gas saturation. We use the term

*gas*to represent a gaseous mixture of uncontaminated soil air and TCE vapors.

In this mathematical model, it was assumed that the gas density was not influenced by the gas pressure because the variation of pressure in time and space would be insignificant in the laboratory column. The model takes into account, however, the influence of the TCE concentration vapors on the gas density.

*C*

_{a}(kg m

^{−3}) is the TCE vapor concentration in the gas mixture,

*M*

_{air}(g mol

^{−1}) and

*M*

_{TCE}(g mol

^{−1}) represent the molecular weights of the uncontaminated soil air and the TCE, respectively, and ρ

_{air}(kg m

^{−3}) is the density of the uncontaminated soil air.

_{a}, is expressed by

_{air}(Pa s) is the dynamic viscosity of the uncontaminated soil air, and

*C*

_{a,sat}(kg m

^{−3}) and μ

_{a,sat}(Pa s) denote the TCE saturation concentration and the dynamic viscosity of the TCE-saturated gas phase, respectively.

*D*

_{a}

^{eff}(m

^{2}s

^{−1}) denotes the effective air diffusion coefficient, α

_{L}(m) is the longitudinal dispersivity,

*v*

_{a}(m s

^{−1}) is the average pore velocity of the gas flow, ρ

_{b}(kg m

^{−3}) is the bulk density of the porous medium, and Kd

_{a}(m

^{3}kg

^{−1}) is the distribution coefficient between the air and the solid phase. The average pore velocity of the gas flow was obtained from the Darcy velocity divided by the gas content in the porous medium.

*D*

_{w}

^{eff}(m

^{2}s

^{−1}) and

*S*

_{w}(dimensionless) denote the effective water diffusion coefficient and water saturation, Kd

_{w}(m

^{3}kg

^{−1}) is the distribution coefficient between the water and the solid phase, and

*H*is the dimensionless Henry's law constant.

To solve the coupled pair of Eq. [2] and [5] for the dry sand or Eq. [2] and [6] for the sand at irreducible water saturation, the software Multiphysics (COMSOL Inc., Burlington, MA) was used. COMSOL Multiphysics is based on a finite-element discretization scheme and can be used to simulate multiphase mass transport in three dimensions. The module used in this study was the earth sciences multiphase module containing a preprogrammed transport equation (esst) and a flow equation without any change of water saturation in space and time (esdl). The preprogrammed flow equation represents only a simplified flow case in which changes in the fluid density are negligible. Therefore, the given preprogrammed equation was modified to match Eq. [2]. Gas transport caused by advection and hydrodynamic dispersion is also not taken into account in the current version of Multiphysics. Because these functionalities already existed for the water phase, the variable gas was simulated as an equivalent water phase and, if necessary, the presence of the irreducible water saturation was taken into account.

Before any test in the laboratory, preliminary numerical studies were performed using simplified equations to predict the duration of the various experiments and to adapt sampling and gas chromatography-electron capture detection analysis methods to the expected TCE concentrations.

## Discussion of the Experimental Results

The diffusion experiments were performed under dry conditions (dry sand fill) and under irreducible water saturation (sand at *S*_{irw}). Three experiments (Exp. A, B, and C) were conducted for dry sand, and only two replicates (Exp. D and E) were done for the *S*_{irw} sand. The main characteristics of the five experiments are summarized in Table 3. The fitted equivalent gas diffusion coefficients were obtained by minimizing the mean squared errors between the measured and simulated TCE vapor concentrations. To take into account the inhomogeneous distribution of experimental data, the method of weighted mean squared errors was used to give the same importance to vapor concentrations measured at a steady state as well as at a transient state of mass transport.

### Initial, Transient, and Boundary Conditions of the Column Experiments

The γ densitometry measurements made under dry conditions showed a high homogeneity of the column sand (Fig. 2 ). Higher deviations from the mean counts value were observed at the outside edges of the sand column. These deviations show the influence of the locking rings of the column modules on the measurements. Neglecting these irrelevant data, the standard deviations of the counts along the column were 3.5% of the measured values (Exp. A), 2.0% (Exp. B), 1.3% (Exp. C), 1.6% (Exp. D), and 2.2% (Exp. E).

Reasonably constant profiles of water saturation were recorded from the γ densitometry measurements conducted under wet sand fill conditions (Fig. 3 ). The mean water saturation in the column was 13.3% (Exp. D) and 13.1% (Exp. E), and the standard deviations were 16% (Exp. D) and 22% (Exp. E). The measured water saturation profiles are stable with time for the two tests conducted at irreducible water saturation (Fig. 3).

The differential pressure variations remained below 1.2 Pa in the last four tests (Fig. 4 ). The inadvertent rise in pressure observed in the upstream cavity at the beginning of Exp. A was a result of TCE volatilization and the initial very small diameter of the needles. Mass transport results of the first experiment were considered only under steady-state conditions. To relieve the overpressure that was initially observed, the column inlets were equipped with needles with a sufficiently large diameter.

Despite the air conditioning in the laboratory where the tests were conducted, the variations in the room temperature were quite high during the overall duration of experiments. In all five experiments, the temperature during the day and night differed by between 1 and 4°C (see Table 3).

### Measured Trichloroethylene Vapor Concentrations

During the experiments, TCE vapor concentrations were measured at the inlet (upstream part) and outlet (downstream part) of the sand column. Figure 5 shows the normalized concentrations as a function of time; *C*_{a,sat} corresponds to the TCE saturation concentration in the gas phase at the average temperature of each experiment.

As shown in Fig. 5, all of the experiments that were conducted on dry sand and at irreducible water saturation were reproducible. Therefore, the temperature variations recorded during each experiment (see Table 3) did not have a significant influence on the normalized vapor concentrations.

In both types of experiments, the TCE saturation concentration was achieved in the upstream cavity only after 5 to 6 h. The TCE vapor concentrations recorded in the downstream cavity, however, were always lower than the saturation concentrations of TCE in the gas phase. The maximum TCE vapor concentrations were nearly 75 and 70% of the saturation concentration for the dry sand and the sand at irreducible water saturation, respectively. Steady**-**state conditions of vapor transport were reached earlier in the dry sand columns (15 h) than in the sand columns at irreducible water saturation (20 h). The experimental results demonstrated that water saturation was having only a moderate effect on the time to breakthrough and on the final TCE concentration in the downstream cavity. The influence of the measured 13% water saturation on the overall mass transfer was analyzed more quantitatively and is discussed below.

### Modeling of Trichloroethylene Breakthrough Curves Considering Diffusion and Mass Exchange between Interfaces

#### Numerical Implementation and Model Parameters

In a first step, the purely diffusive transport of TCE in the water and gas phases in the sand column and the gaseous–solid partitioning or the gaseous–aqueous and aqueous–solid partitioning at equilibrium were considered without taking into account the influence of the vapor density on the gas pressure field. In the upstream and downstream cavities, only diffusive transport of TCE vapors was considered using the free air diffusion coefficient of *D*_{a}^{0} = 8.1 × 10^{−6} m^{2} s^{−1}.

To express the Henry's law coefficient for mass exchange between the water and gas phases, the empirical relation proposed by Munz and Roberts (1987) for chlorinated solvents was used. For partitioning between the aqueous and solid phases, a linear sorption isotherm was chosen (Jury et al., 1990; Mendoza and Frind, 1990; Sleep, 1998; Kim et al, 2001, 2005). The low sorption capacity of TCE on medium-size quartz sand, as reported by Benremita (2002), justified this choice.

_{a}(dimensionless) is the tortuosity of the porous medium toward air.

The Millington–Quirk expression has been used often (Jury et al., 1990; Mendoza and Frind, 1990; Batterman et al., 1995; Conant et al., 1996; Sleep and Sykes, 1989; Sleep, 1998; Bohy et al., 2006) and can be applied to vapors in sandy materials, as demonstrated by Weeks et al. (1982). With the Millington and Quirk model, the effective air diffusion coefficient was set to 2.4 × 10^{−6} m^{2} s^{−1} for the dry sand and 1.5 × 10^{−6} m^{2} s^{−1} for the sand at irreducible water saturation.

_{w}(dimensionless) is the water tortuosity of the porous medium, estimated in the numerical studies to be 0.0062 and 0.0063 for the sintered glass and sand, respectively, and

*D*

_{w}

^{0}(m

^{2}s

^{−1}) represents the free water diffusion coefficient (see Table 2).

The gas or liquid sorption coefficient used for the dry sand fill and the sand at irreducible water saturation, Kd_{a} or Kd_{w}, as well as Henry's law constant are summarized in Tables 1 and 2.

Supposing that there is no density-driven advection in the gas phase (*v*_{a} = 0), Eq. [5] or [6] were solved using the COMSOL Multiphysics code for the given boundary and initial conditions of the column experiments at a constant temperature of 20°C. The grid used in our simulations was a two-dimensional radial symmetrical (*r*–*z*) mesh with triangular elements. The column axis corresponds to *r* = 0, and four subdomains with different physical properties were considered (Fig. 6 ). Those subdomains were: the upstream cavity, the sand, the sintered material, and the downstream cavity. The model was therefore composed of 6000 finite elements. In the simulations, a fixed-concentration boundary was set in the bottom of the upstream cavity at the DNAPL surface at *z* = −0.03 m (see Fig. 6) equal to the TCE saturation concentration *C*_{a,sat} at 20°C. The top of the downstream cavity, located at *z* = 0.148 m, was chosen to correspond to a no-flow boundary. A time step of 500 s was chosen in the numerical simulations.

#### Numerical Results

Figure 7 shows the TCE concentration breakthrough curves calculated in the downstream cavity compared with those measured in Exp. A, B, and C (dry sand) and Exp. D and E (sand at irreducible water saturation). Large differences exist between the measured and the numerically simulated curves. The mean squared errors between the experimental data and the calculated concentrations were 2.1 for the dry sand and 3.9 for the sand at irreducible water saturation (Fig. 7).

To better approach the observed TCE breakthrough curves, the effective air diffusion coefficients were fitted in a second numerical step. Adopting a manual least mean square error procedure, an effective air diffusion coefficient of *D*_{a}^{eff} = 3.4 × 10^{−7} m^{2} s^{−1} was obtained with a mean squared error of 0.34 for the dry sand. For the diffusion experiments with sand at irreducible water saturation, the fitted coefficient obtained with a mean squared error of 1.27 was *D*_{a}^{eff} = 3.0 × 10^{−7} m^{2} s^{−1}. Those fitted gas diffusion coefficients are seven or five times lower than the empirical coefficient given by Millington and Quirk (1961) under the same conditions. Both of the fitted TCE breakthrough curves are shown in Fig. 7. Note that neglecting mass transfer between interfaces will result in lower fitted gas diffusion coefficients, which are about seven or eight times lower than the ones based on Millington and Quirk (1961).

The maximum TCE vapor concentrations measured in the downstream cavity were about 25% to 30% lower than both the saturation concentration in the upstream cavity and the maximum normalized TCE concentration simulated in the two modeling cases. By neglecting the advective part of the vapor mass transport and fitting an effective air diffusion coefficient to come closer to the experimental results, the steady-state conditions of mass transport will not be reproduced and the real effective air diffusion coefficient will always be underestimated.

The second step of numerical simulations of the measured TCE breakthrough curves is discussed below. This second step considered in detail all relevant transport mechanisms in the sand column (i.e., diffusion, advection, and mass exchange between interfaces) and diffusive transport of TCE in the upstream and downstream cavities. The improved modeling approach also took better account of the measured time-variant TCE vapor concentration in the upstream cavity (see Fig. 5).

### Modeling of Trichloroethylene Breakthrough Curves Considering Diffusion, Advection, and Mass Exchange between Interfaces

#### Numerical Implementation and Model Parameters

As in the first modeling approach (see above), a two-dimensional radial symmetrical (*r*–*z*) mesh with triangular elements was used in the simulations. The model was composed of 6000 finite elements, and four subdomains were considered (Fig. 6).

The coupled flow and transport equations were simultaneously simulated in the region of the medium sand and sintered material. Equations [2] and [5] were used for the dry sand; Eq. [2] and [6] were solved for the sand at irreducible water saturation. Coupling was done at each iteration step using Eq. [3] and [4]. These equations strongly link the gas density and the dynamic gas viscosity to the vapor concentration and thus directly influence the gas pressure field and the mean gas flow velocity.

The empirical relation of Millington and Quirk (1961) was used to express the effective gas diffusion coefficient (see Eq. [7]) and the effective water diffusion coefficient (see Eq. [8]). Analogous to the first modeling approach, the empirical relation proposed by Munz and Roberts (1987) was used to express the Henry's law coefficient for mass exchange between the water and gas phases. For partitioning between the aqueous and solid phases, a linear sorption isotherm was used.

*a*

_{a}denotes the tortuosity–connectivity parameter for gas flow in the porous medium.

*n*and the Brooks–Corey parameter λ (Morel-Seytoux et al., 1996):

A wide range of values for the tortuosity–connectivity parameter *a*_{a} is found in the literature. Very often, a value of 0.5 (Parker et al., 1987; Chen et al., 1999; Monlouis-Bonnaire et al., 2004; Emonet, 2003) or 2 (Thomson et al., 1997; Jellali, 2000) has been used. In some studies (e.g., Fischer et al., 1997), the tortuosity–connectivity parameter was adjusted using experimental air permeability measurements. Various values can be found depending on the fitting approach or the soil type, e.g., 1.0 < *a*_{a} < 2.8 for a sand (Tuli and Hopmans, 2004) and 0.87 (sand) < *a*_{a} < 7.7 (loam) (Tuli et al., 2005). In the numerical study, *a*_{a} was chosen equal to 2. In the sensitivity study presented below, how the modification of this parameter influenced the vapor breakthrough curves is discussed.

Using Eq. [9], the relative gas permeability was fixed to 1 for the dry sand and 0.69 for the sand at irreducible water saturation.

Transport of TCE vapor in the two cavities was modeled by solving the gas diffusion equation using an equivalent diffusion coefficient of *D*_{a}^{eq} = 3.5 × 10^{−7} m^{2} s^{−1} to reproduce the measured concentrations at the top of the upstream cavity in a preliminary laboratory study (Cotel, 2008).

At the lower limit of the upstream cavity at *z* = −0.03 m (see Fig. 6), a fixed-concentration boundary was set in the bottom of the upstream cavity at the DNAPL surface equal to the TCE saturation concentration *C*_{a,sat} at 20°C. The upper limit of the downstream cavity at *z* = 0.148 m (Fig. 6) corresponded to a no-flow boundary; an internal concentration boundary condition linked the four model domains (Fig. 6).

To solve the flow equation, a prescribed gas pressure was chosen at the upstream and downstream limits of the porous medium that corresponded to the pressure of the atmosphere at *z* = 0 and 0.096 m (Fig. 6). The initial conditions of the pressure field represented a vertical pressure profile at static equilibrium. A time step of 100 s was chosen.

The model parameters and boundary and initial conditions of the flow and transport models used to simulate the TCE vapor experiments are summarized in Table 4 .

#### Numerical Results

The TCE vapor concentrations calculated in the upstream and downstream cavities are shown in Fig. 8 for the two types of experiments.

Two major conclusions arise from comparing the simulated breakthrough curves with the measured curves. First, at the sampling point of the upstream cavity (*z* = −0.02m), the measured and simulated concentrations are very close. The chosen apparent air diffusion coefficient was adequate to correctly reproduce the overall mass transport in this part of the setup. Second, the stationary normalized TCE vapor concentrations measured in the downstream cavity were underestimated by up to 9% by the numerical model. This result is acceptable because none of the model parameters was fitted at this stage of the simulation. Furthermore, the observed steep increase in concentration is slightly lower in the numerical breakthrough curves. As expected, the obtained numerical results were significantly closer to the observed concentration in this model because this numerical model took into account density-driven advection. For the dry sand, the mean squared error was 0.055. For the sand at irreducible water saturation, the differences were slightly higher and the mean squared error reached 0.106. A slight difference in the shape of the experimental breakthrough curves and the simulated curves was also observed.

The numerical simulations showed that the 13% water saturation of the medium-size sand had a small impact on the overall mass transfer; however, this water saturation resulted in a lower steady-state concentration of TCE vapor in the downstream cavity. An increase in the water saturation resulted in a reduction in both the air tortuosity and the air permeability of the sand. For low water saturations, an increase in water saturation leads to a greater decrease in the effective air diffusion coefficient than the air permeability of the sand using the selected empirical equations for the effective gas diffusion coefficient and relative gas permeability (Eq. [7] and [9]). An increase in water saturation from 0 to 13% resulted in a 37.1% decrease in the effective diffusion coefficient and a 30.9% decrease in the air permeability. This change caused a relative increase of the downward mass flux compared with the diffusive upward mass flux. This increase was a result of density-driven advection and led to a lower equilibrium concentration. The observed difference between the experimental equilibrium concentrations in dry sand and in sand at irreducible water saturation sand was well reproduced by the numerical model.

In a second numerical step, to better approach the observed TCE breakthrough curves, the effective air diffusion coefficients were fitted. For each experiment, the specific effective coefficient was fixed to obtain calculated breakthrough curves as close as possible to the observed curves (Fig. 9 ). For the dry sand, the effective gaseous diffusion coefficients were 3.4 × 10^{−6} m^{2} s^{−1} (Exp. A), 3.7 × 10^{−6} m^{2} s^{−1} (Exp. B), and 2.9 × 10^{−6} m^{2} s^{−1} (Exp. C), with a mean squared error of 0.014. For the column experiments conducted with sand at irreducible water saturation, effective gaseous diffusion coefficients of 1.8 × 10^{−6} m^{2} s^{−1} (Exp. D) and 1.9 × 10^{−6} m^{2} s^{−1} (Exp. E) were found, with a mean squared error of 0.068. The calculated breakthrough curves were in good agreement with the measured data. In the experiments with dry sand, the standard deviation of the fitted effective gaseous diffusion coefficient for three tests was about 11.3%, and in the tests with sand at irreducible water saturation, the standard deviation for two tests was only 5.5%. The fitted effective gaseous diffusion coefficients were close to those found in the scientific literature (Batterman et al., 1996).

The sensitivity study presented below analyzed in detail the influence of selected model parameters on the overall mass transport and helped to determine if a reasonable modification of their values, within the range of plausible values, might explain the differences between the simulated curves and the measured ones. Another aspect to explore numerically concerns the experimentally observed steady-state conditions of transport, where the diffusive vertical TCE vapor flux is compensated by the downward directed mass flux due to density-driven advection. The second part of the sensitivity study thus focused on the impact of variation of both the imposed upstream vapor concentration and the column height on the steady-state vapor concentrations in the downstream cavity.

#### Sensitivity Study

Effect of Selected Model Parameters on the Simulated Breakthrough Curves. In the flow equation (Eq. [2]) and transport equations (Eq. [5] and [6]), the numerical input values of some model parameters can be adapted within a range of reasonable limits. Using the mean gaseous coefficient calibrated for each type of experiment (see above), a series of numerical sensitivity studies was undertaken. First, the harmonic mean value of the intrinsic permeability of the sand and the sintered material, *k**, was varied within the range 5.7 × 10^{−11} to 6.9 × 10^{−11} m^{2}, which corresponds to a relative variation of ±10%. A variation in the dynamic gas viscosity might have the same effect as the intrinsic permeability. It was decided not to change the value of the dynamic gas viscosity because its value is known precisely (Wilke, 1950). Next, the TCE partitioning coefficients between the gas and the sand, Kd_{a}, and between the water and the sand, Kd_{w}, were modified by multiplying and dividing the initial values by 2, within the limits of 6.0 × 10^{−6} and 2.4 × 10^{−5} m^{3} kg^{−1} for Kd_{a} and 1.1 × 10^{−5} and 4.4 × 10^{−5} m^{3} kg^{−1} for Kd_{w}. In a third step, the dimensionless Henry's law constant *H* was varied between 0.30 (Gossett, 1987) and 0.38 (Conant et al., 1996), where the range corresponds to extreme values given in the literature for TCE at a temperature of 20°C and an atmospheric pressure of 0.1 MPa. Finally, the tortuosity–connectivity parameter that appears in the relative gas permeability curve (Eq. [9]) was modified between 0.5 (Parker et al., 1987) and 2.8 (Tuli and Hopmans, 2004). For the sand at irreducible water saturation (*S*_{irw} = 0.13), the chosen range of *a*_{a} corresponds to a variation of relative gas permeabilities between 0.85 and 0.62.

Table 5 summarizes how the four parameters influenced the simulated breakthrough curves. A detailed comparison of the calculated TCE vapor curves with the measured concentration curves is shown in Fig. 10 .

The simulated concentration breakthrough curves were not very sensitive to variations in the dimensionless Henry's law constant and partitioning coefficients but were sensitive to modifications of the intrinsic permeability and the gas tortuosity–connectivity.

In both types of experiments, the increase in permeability enhanced gravity-driven advection and thus decreased the maximum vapor concentration reached under steady-state conditions in the downstream cavity. Decreasing the tortuosity–connectivity parameter *a*_{a} resulted in an increase in the relative gas permeability. This increase favored density-driven advective mass fluxes, which contributed to a decrease in the maximum steady-state concentrations. As the partitioning coefficient increased, the retardation factor also increased, resulting in an overall shift of the curve. Therefore, the time to achieve steady-state transport conditions increased. An increase in the Henry's law coefficient should have the same effect on the breakthrough curve, but the differences between the low and high values of the Henry's law coefficients found in the literature are too small to visualize the resulting differences in the breakthrough curves.

Fitting effective diffusion coefficients is not the only way to reduce the underestimation of the maximum vapor concentration reached under steady-state conditions. One could also have chosen an intrinsic permeability lower than *k** = 5.7 × 10^{−11} m^{2} and, for the sand at irreducible water saturation, a relative gas permeability lower than 0.62. In order to stay within the range of reasonable limits, it is not possible to correctly predict with the numerical model the maximum concentrations that were observed in the experiments. The mean squared error was about 0.023 for the dry sand and 0.077 for the sand at irreducible water saturation. The intrinsic permeability was measured accurately in preliminary column experiments (Cotel, 2008) (Table 1), and further adaptation of this model input parameter beyond its given physical range is not logical.

Small differences in the curve shape remain between the measured and the calculated breakthrough curves. The modification of Henry's law constant and the partitioning coefficients (Kd_{a} and Kd_{w}) did not significantly improve the calibration of the breakthrough curves. Therefore, their initial values were considered to be correctly chosen for the given experimental cases.

*H** (m) denotes the height of the soil column, Ra (dimensionless) is the Rayleigh number (Simmons et al., 2001; Weatherill et al., 2004; Farajzadeh et al., 2007).

To study the influence of both the Rayleigh number and the upstream boundary concentration on the TCE vapor concentrations in the column outlet, the numerical results of the second step described above were used as references. Two series of numerical simulations were conducted with COMSOL Multiphysics. In the first analysis, the same initial model parameters were used as before, and the imposed upstream concentration (*C*_{a,up}) was varied. Figure 11 shows the numerical results of seven studies for both experimental conditions (dry sand and sand at *S*_{irw}). The upstream vapor concentrations chosen were 80, 60, 40, 20, 10, 5, and 1% of the TCE saturation concentration in the gas phase at 20°C (*C*_{a,sat} = 0.419 g L^{−1}). As illustrated in Fig. 11c, the lower the vapor concentration at the upstream boundary, the lower the influence of the gas density, ρ_{a}, on the normalized maximum downstream vapor concentration. Assuming that density-driven advection is negligible when the downstream vapor concentration achieves 95% of the prescribed upstream vapor concentration, it can be noted that this may be the case for a relative gas density of 1.044 and 1.034 for the dry sand and sand at *S*_{irw}, respectively, based on the results shown in Fig. 11c. When Δρ = (ρ_{a} − ρ_{air}) tends to zero, the limit case where *C*_{a,down} ∼ *C*_{a,up} is obtained. This corresponds to a Rayleigh number of zero (see Eq. [11]).

In the second series of numerical simulations, only the height of the soil column was increased. In this case, the vapor concentrations at the upstream boundary were kept constant and equal to the TCE saturation concentration in the gas phase at 20°C (*C*_{a,sat} = 0.419 g L^{−1}). The heights of the simulated soil columns were about 1.6, 2.1, 4.2, 6.3, 8.4, and 10.4 times the real column height of 0.096 m. The influence of the column height on the steady-state vapor concentrations at the column outlet is illustrated in Fig. 12 . The higher the column, the lower the maximum vapor concentration in the column outlet. This is not really surprising. The quantified relationship, however, is not linear: for example, in the case of the dry sand, doubling the height of the column led to a normalized maximum vapor concentration of about 0.58, whereas in the case of a four times higher column, the maximum ratio achieved 0.40. This shows that with increasing Rayleigh numbers, mass flux due to density-driven advection becomes dominant and thus significantly reduces vapor diffusion in the vertical direction.

## Conclusions

The transport of TCE vapors in low water-saturated porous media was studied using soil column experiments with a medium-size sand at two different water saturations, completely dry and irreducible water saturation.

The measured vapor breakthrough curves were reproducible for both water saturations tested. The maximum TCE vapor concentrations were 7% higher in the dry column. These concentrations were only 75% of the prescribed inlet vapor concentration, indicating the strong influence of the density-driven advective mass transport on the vertical diffusion of TCE. The experimental results highlight that water saturation has only a moderate effect on the time to breakthrough and the final TCE concentration in the downstream cavity. In a field scenario, however, the conclusion might not be the same for two reasons. First, the effective gas diffusion coefficients might be much smaller in the field due to high water saturation in the neighborhood of the capillary fringe. Second, mass transfer between gas and water phases, especially as the amount of water compared with the gas volume containing TCE vapor increases, cannot be described by a gaseous–liquid partitioning at equilibrium and needs to be expressed by a kinetic transfer term.

Without taking into account the density of the vapors in the soil air, the derived apparent TCE diffusion coefficient was lower than the coefficients given in the literature by a factor of seven to eight in both cases when vapor diffusion was the only considered transport mechanism. For the given experimental conditions, however, diffusion of TCE in the residual water saturation combined with mass transfer between solid, gaseous, and water phases had an influence on the transient part of the measured vapor breakthrough curves and thus cannot be neglected. It is worthwhile to note that by considering all the mechanisms mentioned above, the derived apparent TCE was still lower than the coefficients given in the literature by a factor of five to seven. This clearly highlights the significant influence of the TCE vapor density under the given experimental conditions.

The results of the numerical study underscore that an exact representation of the dimensionless Henry's law constant is not essential to accurately reproduce the transport of TCE under conditions similar to those in the presented experiments. Conversely, the intrinsic permeability of the porous medium and the gas tortuosity–connectivity parameter were found to be important parameters that influenced the predicted vapor concentrations. Density-driven mass flux may play a significant role, especially in natural soils with similar high permeabilities, and predictions of the fate of vapor plumes must be based on an accurate evaluation of the permeability field within the vadose zone.

In natural soils, strong density gradients may exist within the vapor plume and at its fringe, resulting in a significant advective transport of DNAPL vapors. Using vapor concentration gradients of the plume to evaluate mass fluxes based on Fick's law, without taking into account the effect of density-driven advection, will therefore lead to an overestimation of the upward-directed vertical total mass flux of the vapors. From the sensitivity study, interesting lessons can be learned. For the given 22-cm-long medium-sand column at irreducible water saturation, downstream TCE vapor concentrations achieved nearly 95% of the upstream vapor concentrations, maintained constant at about 10% of the TCE saturation concentration. In this case, a maximum relative gas density of about 1.034 was quantified at the upstream boundary of the medium-size sand and vapor migration due to density gradients could be neglected. An increase in height of the soil column, however, corresponding to an increase in the Rayleigh number, resulted in a significant decrease in the overall vertical mass flux as density-driven advection became dominant compared with diffusion. Combining both results, it can be concluded that even for gas samples taken in the field, with concentrations <10% of the TCE saturation concentration, density-driven advection may represent a significant transport mechanism for vapors in the unsaturated zone.

## Appendix A

### Steady-State Trichloroethylene Vapor Flux

_{diff,z}(kg s

^{−1}), in the soil column is equilibrated by the downward-directed TCE mass flux due to gravity-driven advection, ṁ

_{diff,z}(kg s

^{−1}):

*A*(m

^{2}) is the cross-sectional area of the soil column.

*C*

_{a,down}) and the upstream boundary (

*C*

_{a,up}), the concentration gradient can be written as

*H**(m) denotes the height of the soil column.

*z*direction was obtained from the Darcy velocity divided by the gas content of the porous medium:

*h*=

*p*

_{a}/(ρ

_{air}

*g*) +

*z*introduced by Lusczynski (1960), the Darcy velocity in the

*z*direction may be expressed by

*h*. Given that ∂

*h*/∂

*z*is zero, Eq. [5A] simplifies thus to

Financial support for this research from the Agence de l'Environnement et de la Maîtrise de l'Energie (ADEME) and the Commissariat à l'Energie Atomique (CEA) is gratefully acknowledged.