CHAPTER 2 METHEMATICAL MODEL
2.1 Mathematical model: Two-component case
2.1.1 Below the evaporation front
In this region, NAPL phase exists and the saturation of NAPL phase is greater than zero, i.e., SR >0. Equation (8) can then be expressed as
with representing mole fractions of organic compounds in the NAPL phase, i.e., . For a one-dimensional system, Equation (6) becomes:
ui based on Equation (11), the mass-conservation equation for each component can be written
P Gi R i ρ φ/C σ =
as: NAPL saturation S0 and mole fraction u0 of component one as:
)
Substituting Equations (10) and (14) into Equation (12) yields
( )
2.1.2 Above the evaporation frontIn this region the NAPL phase fully evaporates. i.e., SR =0. For a two-component VOC, Equation (11) can be written for each component as follows:
2
where u1 and u2 represent normalized concentrations above the evaporation front,
2.1.3 At the evaporation front
Assume that the evaporation of NAPL phases occurs instantaneously and the sum of mole fraction equals one, i.e., u1+u2 = 1. Combining Equations (12) and (13), the equation for mass conservation at the front can be expressed as
t
Taking the limits of and , Equation (18) describing the front can then be written as
where the superscripts t+ and t– of the mole fraction in Equation (20) denote the mole fraction at the time slightly after and before the volatilization, respectively, and z+ and z– represent the mole fraction at the region slightly below and above the front, respectively. Equation (20) can be further simplified as
(
t t) (
t t)
R2.1.4 Boundary and initial conditions
For simulating the ambient environment near the surface, the gas phase VOC is
considered to diffuse across a stagnant air boundary layer with thickness d and the gas phase concentration is assumed equal to zero at the top of boundary layer. The flux diffusing to the atmosphere can then be expressed as (Jury et al., 1983)
G T G G
R hC z
D C =
∂
∂ , at z=0 (22)
where . If the thickness of air boundary layer near the surface is very small and negligible, one may assume that d = 0. Equation (22) can then be expressed as
d D h= Gair /
=0
CT , at z=0 (23)
For a multicomponent VOC, the upper boundary conditions can be written from Equation (23)
as
( )
0, 2( )
0, 01 t =u t =
u (24)
The conditions of the mole fractions at the lower boundary are denoted as:
( )
L t uou1 , = , u2
( )
L,t = 1−uo (25) where L is the depth of the lower boundary.The initial NAPL saturation and the mole fractions of each component are assumed
spatially uniform; they are
( )
z,0 S0SR = (26)
( )
z uou1 ,0 = , u2
( )
z,0 =1−uo (27)The evaporation front is initially located at the land surface and thus denoted as
( )
0 =0s (28)
2.2 The numerical method in solving the model
This section presents the numerical method used to solve the two-component model.
2.2.1 Finite difference approximation
The equations of describing VOC transport in the three regions are solved separately by the finite difference method. An interpolative moving grid approach (Javierre et al. 2006) is adopted to handle the moving boundary problem. The total number of nodes within the problem domain is equal to N and r is the nodal number assigned at the evaporation front beginning from the land surface. Accordingly, the number of grids from the land surface to the front is r-1, the number of grids below the front is N-r, and the initial grid size dz is equal to L/(N-1). The grid sizes above and below the front defined as dzr and dzN-r, respectively, need to be re-estimated after each move of the front. To avoid introducing large truncation error, the grid sizes dzr and dzN-r should be close to dz. If the front move to a location between the nodal numbers initially assigned as j and j+1, the new grid sizes of dzr and dzN-r
are then estimated by s/(r-1) and (L-s)/(N-r), respectively, where r = j+1. The backward difference relative to time and central difference relative to space are used to approximate the VOC transport equations. Therefore, the difference equation for the mole fraction in the
region below the front obtained from Equation (16) is
where n is the number of time step, and dt is the time interval. The difference equations for the mole fractions of the two components in the region above the front obtained from Equations (17) and (18) are, respectively:
2
Finally, the difference equations for the mole fractions of the two components at the front also derived from Equations (17) and (18) with different grid sizes above and below the front are, respectively, expressed as
r
Substituting Equation (10) into Equation (21), the difference equation for the front location can then be obtained as
( ) (
n)
R2.2.2 The solution procedure of the model
Because the location of the evaporation front s(t) is unknown, a trial and error procedure is therefore taken to find the front location. The steps of the procedure are listed below and the related flowchart is shown in Figure 2:
1. Give the initial front location (Equation (28)) and nodal values of mole fraction based on the boundary conditions (Equations (24) and (25)) and initial conditions (Equations (26) and (27)).
2. Guess front location after the start of evaporation.
3. Determine the nodal number for the front and the grid sizes based on the front location.
4. Solve the solutions for the region below the front (Equation (29)), for the region above the front (Equations (30) and (31)), and at the front (Equations (32) and (33)).
5. Compute the front location based on Equation (34).
6. Proceed to next time step if the location of the front converges, i.e., the difference between two succeeding estimations of the front location is less than a very small value, e.g., 10-10 m;
otherwise, go back to step 3.
2.3 Simplified model: Single-component case
This section presents a single-component model with the governing equation simplified from the two-component model. The domain of the single-component model is also divided into three regions based on the front location. In a one-dimensional homogeneous and unsaturated soil system, Equation (6) representing the VOC transport can be written as
2 0 2.3.1 Below the evaporation front
In this region, the NAPL has not vaporized to gas phase yet. The VOC concentrations in each phase are the initial saturated concentrations and the total concentration can be expressed as
is the saturated VOC concentrations for pure component in adsorbed phases, CT0
is the initial total concentration, and θG0, θL0, and θR0
are the initial volumetric contents, respectively.
2.3.2 Above the evaporation front
In this region, NAPL phase of the VOC completely vaporizes to gas phase; therefore, the VOC presents only in gas, liquid, and adsorbed phases, i.e., CT = CGθG+ CLθL+CSρb. Jury et al. (1983) used a ratio R to represent each phase in relation to the total concentration.
Accordingly, one may introduce the following equation based on Equations (1) and (3):
With Equation (37), Equation (35) can be written as
2 0 where denotes as effective diffusion coefficient. Equation (38) describes the VOC transport in gas and liquid phases between the land surface and evaporation front. In reality, the volatilization of NAPL occurs right at the evaporation front which will be discussed in the next section.
G G
E D R
D = /
2.3.3 At the evaporation front
The evaporation front moves downward while the NAPL vaporizes to the gas phase.
Assume that the transformation of the contents between these two phases occurs instantaneously and the liquid volumetric content is unchanged, i.e., θL = θL0
. In addition, the VOC concentrations in each phase are still saturated,i.e., CG = CGP
, CL = CLP
, and CS = CSP
. The total concentration at the front can therefore be expressed as
0 T
T C
C = (39)
The gas, liquid, and NAPL phases remain in the unsaturated soil pores, the sum of θG, θL, and
θR equals soil porosity φ, i.e., θG+θL+θR =φ. Assumes φ does not change with time,
the transformation of volumetric content with respect to time among each phase is conserved.
Thus,
following relationship from Equation (1)
( )
the front can then be written as( )
where the superscripts + and - denote the volumetric content at the time slightly after and before the volatilization, respectively, and the concentrations at the region slightly below and above the front, respectively. Consider that the volatilization occurs instantaneously, therefore equals and equals zero after evaporation. The VOC
concentrations in each phase are the initial saturated concentrations below the front and the concentration gradient of gas phase below the front is naturally equal to zero, i.e.,
+
θR θR0 θR−
=0
∂
∂CG+ z . Since the liquid density is higher than the gas phase concentration about three
orders (Falta et al., 1989, Tables 1 and 4), the term related to CGP
on the left-hand side of Equation (42) is thus negligible. Accordingly, Equation (43) representing the front z = s(t) can be expressed as
z
2.3.4 The analytical solution of single-component model
Consider that the VOC is saturated or in equilibrium state in different phases and uniformly distributed in the unsaturated soil initially. The mathematical model describing the single-component VOC transport in the soil consists of Equation (38) as the governing equation, Equations (39) and (43) as the lower boundary conditions, and Equation (23) as the upper boundary condition.
Based on Boltzmann’s transformation, a new variable is defined as ξ = z 2 DEt . Equation (38) can then be transformed to an ordinary differential equation as
0
The solution of Equation (44) can be obtained as (Carslaw and Jaeger, 1959)
( )
B erfA
CT(ξ)= ⋅ ξ + (45)
where erf(ξ) is the error function and A and B are unknown coefficients. Substituting Equation (45) into Equation (23), the result for the concentration distribution is
( )
⎟⎟Substituting Equation (39) into Equation (46), the evaporation front s(t) and coefficient A can then be obtained, respectively, as
t
where α is an unknown constant depending upon the soil parameters and contaminant characteristics.
The time of vanish of NAPL can be solved by Newton’s method (Yeh, 1987) from Equation (47) when the front reaches a target location below the land surface designated by the environmental or legal requirement. In addition, the moving speed of the evaporation front can also be obtained after taking the derivative of Equation (47) with respect to time and the result is
Uf t 2
= α (49)
Substituting Equations (46), (47) and (48) into Equation (43) yields
⎟⎟
The unknown constant α can then be easily determined from Equation (50) by Newton’s method. Note that the normalized total concentration is defined as CT(z,t)/CT0
, representing the mole fraction in the single-component model.
CHAPTER 3 RESULTS AND DISCUSSION
Leaks of petroleum fuels from the underground storage tanks are common problems for soil contamination. The petroleum spills are often associated with aromatic hydrocarbons such as benzene, toluene, ethyl benzene, and various xylene isomers (BTEX). In this section, the hydrocarbons of benzene and toluene are chosen to simulate their transport and mole fraction distributions in unsaturated soils using the two-component model. In the past, Carbon tetrachloride was commonly used as coolant in industry or produced as the fire
extinguishers. Carbon tetrachloride is highly toxic; a small amount of this chemical residing in the soil may pose severe environmental problems. The carbon tetrachloride in the
unsaturated soil is considered as a target VOC and analyzed using the single-component model.
Six cases are considered to address the issues in regard to the evaporation rate,
evaporation front movement, mole fraction, and concentration distributions of VOC for the present models. Case 1 is to compare the mole fractions of toluene at various evaporation times predicted by single-component and two-component models. Case 2 examines the effect of initial mole fraction on the evaporation and the changes of the mole fraction distributions of benzene and toluene. Case 3 investigates the effect of soil porosity on vaporization of carbon tetrachloride from NAPL phase to gas phase while case 4 addresses the issue of carbon tetrachloride transport based on the model with and without considering
the presence of NAPL phase. Case 5 studies the migrations of evaporation front for different contaminants, namely carbon tetrachloride and toluene. Case 6 discusses the effect of
effective diffusion coefficient of carbon tetrachloride on the moving speed of evaporation front. The values of the soil chemical properties are listed in Table 1 and the properties of benzene, toluene, and carbon tetrachloride are given in Table 2 for these six cases. Note that the depth of the lower boundary L is chosen as 5 m, the total number of nodes N is 10000, and the time interval dt is 0.1 sec in the case study when adopting the finite different
approximation for the two-component model.
3.1 Case 1: Different evaporation times in two models
This case uses the same assumptions for both single-component and two-component models and considers that toluene is the only VOC found in the soil, i.e., u0 = 1. Figure 3 shows the mole fraction distributions of toluene versus depth predicted by the
single-component and two-component models at various evaporation times. The dashed line denotes the solution of single-component model while the solid line represents the results predicted by the two-component model. Moreover, the symbols of rhombus, triangle, and circle represent the mole fractions at times 1, 10, and 100 day, respectively. This figure shows the front locations at various evaporation times and at the front the mole fraction equals its initial value for the single-component VOC. The figure also shows that the curves
predicted by both models are fairly close, implying that the results predicted by the
two-component model with the present numerical approach match well with those estimated based on the analytical solution of the single-component model. The moving speeds of the front Uf estimated by equation (49) are 8.296×10-2, 2.624×10-2, and 8.296×10-3 m/day at times 1, 10, and 100 day, respectively, indicating that the moving speed decreases rapidly at early time and then slowly as time increases
3.2 Case 2: Initial mole fraction
In this case, benzene is considered to be component one and toluene is component two in the two-component model. Figures 4(a) - 4(c) show the mole fraction distributions of
benzene and toluene versus depth when the initial mole fractions of component one are 0.2, 0.5, and 0.8, respectively, at 100 day. The evaporation front of the NAPL with u0 = 0.2, 0.5, and 0.8 reaches 0.860, 0.931, and 1.002 m below the surface, respectively. In addition, at the front u1 = 0.123 and u2 = 0.877 when u0 = 0.2, u1 = 0.310 and u2 = 0.690 when u0 = 0.5, and u1 = 0.498 and u2 = 0.502 when u0 = 0.8. The figures show that the depth of the front increases with the initial mole fraction of benzene, representing the moving speed of the front depends on the initial mole fraction. The mole fraction of benzene increases with depth until reaching u1 = u0; on the other hand, the mole fraction of toluene increases above the front but decreases below the front until reaching u2 =1 u− 0. These results indicate that at the front, the mole fraction of benzene decreases as time increases while that of toluene increases with time. Moreover, the mole fractions of both components are not equal to their initial values at
the front. In fact, both components reach their initial values occurred at certain distances below the front. Such a phenomenon can be attributed to the fact that benzene has higher evaporation efficiency than toluene. Therefore, the mole fractions of benzene and toluene change with depth, although the NAPL below the front does not evaporate. Figure 5 shows the curve of NAPL phase saturation, calculated from Equation (14), versus mole fraction of benzene. The figure demonstrates that SR equals 0.004 when u0 = 0 and 0.015 when u0 = 1.
This result indicates that SR changes with mole fraction and varies slightly below the front.
3.3 Case 3: Soil porosity
This case examines the effect of soil porosity on the concentration distribution of carbon tetrachloride using the single-component model. The total evaporation time is considered to
be 100day. Note that the saturation S of each phase is constant and the volumetric contents θ of gas, liquid, and NAPL change in equal proportion with the soil porosity. Figure 6
shows the predicted normalized total concentration versus depth at 100 day for the porosities of 0.1, 0.3, and 0.5. It is apparent from Figure 4 that the vaporization increases moderately with soil porosity. In addition, the evaporation front of the NAPL with φ = 0.1, 0.3, and 0.5 are at 0.831, 1.432, and 1.860 m below the land surface, respectively, and the moving speeds of the front are 4.154×10-3, 7.158×10-3, and 9.301×10-3 m/day, respectively, indicating that the moving speed of the front increases with soil porosity although different amounts of NAPL exist in the soil.
3.4 Case 4: Absence of the NAPL phase
In this case, the effects of presence and absence of NAPL on the distribution of carbon tetrachloride in the soil are compared and studied. Jury et al. (1983) considered the scenario that there were three phases of VOC presented in the soil with neglecting the NAPL phase.
Their analytical model included the mechanisms of diffusion, soil water advection, and first-order decay. They used a diffusive flux as the upper boundary condition at the land surface and zero total concentration at infinite depth as the lower boundary condition. Jury et al.’s model is simplified by neglecting the water phase advection and decay and thus called simplified Jury et al.’s model hereafter.
Figure 7 shows the curves of normalized total concentration of carbon tetrachloride versus depth predicted by two different VOC transport models at t = 100 day. The solid line with triangle symbol represents the normalized total concentration distribution predicted by the present single-component model while the solid line with circle predicted by simplified Jury et al.’s model. Figure 7 indicates that the normalized total concentration of VOC predicted by the present model is significantly higher than that of simplified Jury et al.’s model. Although NAPL occupies only one-percent of volume in the soil pores, the NAPL however affects the total concentration distribution and transport capability greatly. Figure 8 exhibits the predicted distribution of normalized total concentration versus time for the
present model and simplified Jury et al.’s model at the depth of 1 m. The figure shows that
the total concentration of carbon tetrachloride decreases quickly after 2 day predicted by simplified Jury et al.’s model and after 37 day by the single-component model indicating that the presence of NAPL has significant impact on the VOC transport.
3.5 Case 5: Different chemicals
In this case, both carbon tetrachloride and toluene are considered to reside in the unsaturated soil. Table 2 shows that toluene has less molecular weight and liquid density and lower saturated vapor pressure and Henry’s law constant than those of carbon
tetrachloride. Figure 9 shows the curves of normalized total concentration versus depth for carbon tetrachloride and toluene at 100 day predicted by the single-component model for each chemical. The solid lines with circle and triangle represent the concentration distributions of carbon tetrachloride and toluene, respectively. The vaporization of toluene is significantly lower than that of carbon tetrachloride; the depths of the evaporation front of carbon
tetrachloride and toluene equal 1.659 m and 0.813 m, respectively, at 100 day. Figure 10 shows that the depths of the front versus evaporation time for both chemicals. This figure indicates that the fronts at 50 day and 100 day reach the depths of 1.173 m and 1.659 m, respectively, for carbon tetrachloride and 0.575 m and 0.813 m, respectively, for toluene.
Obviously, the migration of the front of carbon tetrachloride is significantly faster than that of toluene due to lower vapor pressure and Henry’s Law constant value of toluene. Such a phenomenon results in lower vapor pressure gradient and proportion of concentration in gas
and liquid phases. In addition, the diffusion coefficients of gas and liquid phases of both chemicals differ by four orders of magnitude. Thus the diffusion of liquid phase is significantly lower than that of gas phase.
3.6 Case 6: Effect of effective diffusion coefficient on moving speed of evaporation front
Equation (49) shows the moving speed of evaporation front Uf which in fact represents the evaporation rate of NAPL in soil. This case investigates the effect of the soil chemical properties on Uf based on the single-component model. The DE is a function of soil porosity, volumetric content of each phase, air diffusion coefficient, soil bulk density, Henry’s Law constant, and organic carbon fraction. Obviously, different soil chemical properties will affect the value of DE. Figure 11 shows the curves of the depth and moving speed of the front versus evaporation time for DE = 10-10, 10-9, 10-8, and 10-7 m2/s. The depth of the front increases with time and DE greatly while the moving speed is maximal when VOC begins to evaporate and then decreases with increasing time for a fixed value of DE.