Prediction of Vertical DNAPL Vapour Fluxes in Soils Using Quasi-Analytical Approaches: Bias Related to Density-Driven and Pressure-Gradient-Induced Advection

Abstract

This study focuses on a detailed analysis of the errors introduced by two quasi-analytical approaches based on either Fick’s first law or a combination of Fick’s and Darcy’s laws to evaluate the vapour flux of chlorinated solvents from a source zone located in the unsaturated zone towards the atmosphere. A coupled one-dimensional numerical flow and transport model was developed and applied to three case studies characterised by different water content profiles in the vadose zone and under different levels of maximum dense nonaqueous-phase liquid vapour concentrations and vapour pressure conditions of the source zone. The steady-state concentration and pressure profiles obtained were then used in the two quasi-analytical approaches to estimate the flux towards the atmosphere. When mass fluxes due to density-driven advection become dominant and the vertical advective mass fluxes are increased due to strong pressure gradients in the soil air, the error was observed to increase when using the pure diffusion approach in the quantification of the surface flux calculated by the numerical model with increasing dimensionless Rayleigh numbers. Without taking into account the advective transport in the approach, the relative error calculated with only Fick’s law overestimates the real vapour flux when density-driven advection is dominant and underestimates it when pressure-gradient-driven advection dominates. The more advanced advective–diffusive quasi-analytical approach fits reasonably well with the numerically obtained mass fluxes except near soil layer discontinuities, where the evaluation of both the concentration gradient and pressure gradient in the porous media as well as the determination of the average effective diffusion coefficients are rendered more difficult.

Appendices

Appendix A

The flow of soil air is given by the mass conservation equation, which contains mass flux divergence and variation of storage terms as a function of time:

$$ - \nabla \left( {\rho v} \right) = \frac{{\partial \left( {n\rho } \right)}}{{\partial t}} $$

(27)

where ρ (in kilogrammes per cubic metre) is the vapour density, v (in metres per second) is the Darcy’s velocity and n (−) the porosity of the porous media.

The storage term can be explained as follows:

$$ \frac{{\partial \left( {n\rho } \right)}}{{\partial t}} = n\frac{{\partial \rho }}{{\partial t}} + \rho \frac{{\partial n}}{{\partial t}} $$

(28)

Many assumptions are considered to express the flow equation. The first one is to assume that if the compressibility of the medium $\frac{{\partial n}}{{\partial t}}$ is negligible compared with the gas compressibility$\frac{{\partial \rho }}{{\partial t}}$, the storage term will be simplified and defined as follows:

$$ \frac{{\partial (n\rho )}}{{\partial t}} = n\frac{{\partial \rho }}{{\partial t}} $$

(29)

Introducing the ideal gas state equation $\rho \left( {p,T} \right) = \frac{{pM}}{{RT}}$ into Eq. (29), the term of storage variation becomes equal to

$$ n\frac{\partial }{{\partial t}}\left( \rho \right) = \frac{{nM}}{{RT}}\frac{{\partial p}}{{\partial t}} $$

(30)

where M (in kilogrammes per mole) is the molecular weight of the gas, p (in Pascals) is the pressure, R (8.314 Pa m³ mol⁻¹ K⁻¹) is the universal ideal gas constant and T (in Kelvins) is the absolute temperature.

In the divergence term, the velocity of the TCE vapour v is given by Darcy’s law:

$$ v = - \frac{k}{\mu}\nabla \left( {\rho gz + p} \right) $$

(31)

where k (in square metres) is the gas permeability, g (in metres per square second) is the gravity acceleration, μ (in kilogrammes per metre per second) is the dynamic viscosity of the fluid and z (in metres) is the elevation. The left-hand side term of Eq. (30) can thus be rewritten as:

$$ - \nabla \left( {\rho v} \right) = - \nabla \left( { - \frac{{k\rho }}{\mu}\nabla \left( {\rho gz + p} \right)} \right) $$

(32)

$$ - \nabla \left( {\rho v} \right) = \nabla \left( {\frac{{k\rho }}{\mu }gz\nabla \left( \rho \right) + g\rho \nabla (z) + \nabla (p)} \right) $$

(33)

$$ = \nabla \left( {\frac{{k\rho }}{\mu }gz\nabla \rho } \right) + \nabla \left( {\frac{{k\rho }}{\mu }g\rho \nabla z} \right) + \nabla \left( {\frac{{k\rho }}{\mu}\nabla p} \right) $$

(34)

Assuming that variation of the vapour density ρ (p, T) with elevation z is negligible compared with the variation with pressure (p) or temperature (T), one can neglect the first term of Eq. (34). Using the ideal gas state equation, Eq. (30) takes the following developed form:

$$ - \nabla \left( {\rho v} \right) = \nabla \left( {\frac{{k{M^2}{p^2}}}{{\mu {R^2}{T^2}}}g} \right) + \nabla \left( {\frac{{kMp}}{{\mu RT}}\nabla p} \right) $$

(35)

M, g, R, T being constants, Eq. (35) can be rewritten:

$$ - \nabla \left( {\rho v} \right) = \frac{{{M^2}g}}{{\mu {R^2}{T^2}}}\nabla \left( {k{p^2}} \right) + \frac{M}{{\mu RT}}\nabla \left( {kp\nabla p} \right) $$

(36)

The term p∇p may be replaced by

$$ \left( {p\frac{{\partial p}}{{\partial z}}} \right) = \frac{1}{2}\left( {p\frac{{\partial p}}{{\partial z}} + p\frac{{\partial p}}{{\partial z}}} \right) = \frac{1}{2}\frac{{\partial \left( {p.p} \right)}}{{\partial z}} = \frac{1}{2}\frac{{\partial {p^2}}}{{\partial z}} $$

(37)

Using Eq. (37) in Eq. (36) leads to:

$$ - \nabla \left( {\rho v} \right) = \frac{{{M^2}g}}{{\mu {R^2}{T^2}}}\frac{\partial }{{\partial z}}\left( {k{p^2}} \right) + \frac{M}{{2\mu RT}}\frac{\partial }{{\partial z}}\left( {k\frac{{\partial {p^2}}}{{\partial z}}} \right) $$

(38)

Based on Eqs. (30) and (38), one can express the mass conservation equation (see Eq. 27) of 1D vertical gas flow in a porous medium based on pressure p:

$$ \frac{{nM}}{{RT}}\frac{{\partial p}}{{\partial t}} = \frac{{{M^2}g}}{{\mu {R^2}{T^2}}}\frac{\partial }{{\partial z}}\left( {k{p^2}} \right) + \frac{M}{{2\mu RT}}\frac{\partial }{{\partial z}}\left( {k\frac{{\partial {p^2}}}{{\partial z}}} \right) $$

(39)

To linearise Eq. (39), the term p ² is replaced by p using the following expression:

$$ \frac{{\partial {p^2}}}{{\partial t}} = \frac{{\partial \left( {p.p} \right)}}{{\partial t}} = p.\frac{{\partial p}}{{\partial t}} + p\frac{{\partial p}}{{\partial t}} = 2p.\frac{{\partial p}}{{\partial t}} $$

(40)

Replacing $\frac{{\partial p}}{{\partial t}}$ by $\frac{{1\,\partial {p^2}}}{{2\;p\;\partial t}}$ in Eq. (39), one obtains

$$ \frac{{\partial {p^2}}}{{\partial t}} = 2\frac{{Mp}}{{n\mu RT}}g\frac{\partial }{{\partial z}}\left( {k{p^2}} \right) + \frac{p}{{n\mu }}\frac{\partial }{{\partial z}}\left( {k\frac{\partial }{{\partial z}}{p^2}} \right) $$

(41)

As soil gas may occupy only a part of the given pore volume, one must replace porosity n by the gas content θ _g , leading to the following equation:

$$ \frac{{\partial {p^2}}}{{\partial t}} = \frac{{{p}}}{{\mu {\theta_{\mathrm{g}}}}}\frac{\partial }{{\partial z}}\left( {k\frac{{\partial {p^2}}}{{\partial z}}} \right) + 2\frac{{Mg{p}}}{{\mu RT{\theta_{\mathrm{g}}}}}\frac{\partial }{{\partial z}}\left( {k{p^2}} \right) $$

(42)

To simplify the notation of Eq. (42), one sets P = p ² and finally obtains:

$$ \frac{{\partial P}}{{\partial t}} = \beta \frac{\partial }{{\partial z}}\left( {k\frac{{\partial P}}{{\partial z}}} \right) + 2\gamma \frac{\partial }{{\partial z}}\left( {kP} \right) $$

(43)

where $\beta = \frac{{{P_0}}}{{\mu {\theta_{\mathrm{g}}}}},\gamma = \frac{{Mg{P_0}}}{{\mu RT{\theta_{\mathrm{g}}}}}$and p ₀ represents the initial pressure in the soil gas.

Appendix B

The integrals of the linearised flow equation (Eq. 8) are defined as follows:

$$ \int\limits_{{k - 1/2}}^{{k + 1/2}} {\frac{{\partial P}}{{\partial t}}} {\mathrm{d}}z = \int\limits_{{k - 1/2}}^{{k + 1/2}} {\beta \frac{\partial }{{\partial z}}} \left( {k\frac{{\partial P}}{{\partial z}}} \right){\mathrm{d}}z + \int\limits_{{k - 1/2}}^{{k + 1/2}} {2\gamma \frac{\partial }{{\partial z}}} \left( {kP} \right){\mathrm{d}}z $$

(44)

Using an implicit time scheme (discretisation at time t + 1), the integrals for cell k can be expressed as:

$$ P_{{k - 1}}^{{t + 1}}{\mathrm{GP}}(k) + P_k^{{t + 1}}{\mathrm{EP}}(k) + P_{{k + 1}}^{{t + 1}}FP(k) = P_k^t $$

(45)

where

$$ {\mathrm{EP}}(k) = \left[ {1 + {\beta^k}\frac{{\varDelta t}}{{\varDelta {z^2}}}{k^{{k - 1/2}}} + {\beta^k}\frac{{\varDelta t}}{{\varDelta {z^2}}}{k^{{k + 1/2}}} - {\gamma^k}\frac{{\varDelta t}}{{\varDelta z}}{k^{{k + 1/2}}} + {\gamma^k}\frac{{\varDelta t}}{{\varDelta z}}{k^{{k - 1/2}}}} \right] $$

(46)

$$ {\mathrm{FP}}(k) = \left[ { - {\beta^k}\frac{{\varDelta t}}{{\varDelta {z^2}}}{k^{{k + 1/2}}} - {\gamma^k}\frac{{\varDelta t}}{{\varDelta z}}{k^{{k + 1/2}}}} \right] $$

(47)

$$ {\mathrm{GP}}(k) = \left[ { - {\beta^k}\frac{{\varDelta t}}{{\varDelta {z^2}}}{k^{{k - 1/2}}} + {\gamma^k}\frac{{\varDelta t}}{{\varDelta z}}{k^{{k - 1/2}}}} \right] $$

(48)

GP(k) is the factor accounting for the mass flux in the upstream cell k − 1, EP(k) is a factor representing the mass flux in cell k and FP(k) is mass flux in the downstream cell k + 1.

Eq. (45) is given in matrix form as follows:

where

$$ {\text E}P(1) = \left[ {1 + 2{\beta^1}\frac{{\varDelta \tau }}{{\varDelta {\zeta^2}}}{\kappa^1} + {\beta^1}\frac{{\varDelta \tau }}{{\varDelta {\zeta^2}}}{\kappa^2} - {\gamma^1}\frac{{\varDelta \tau }}{{\varDelta \zeta }}{\kappa^2} + {\gamma^1}\frac{{\varDelta \tau }}{{\varDelta \zeta }}{\kappa^1}} \right] $$

(49)

$$ FP(1) = \left[ { - {\beta^1}\frac{{\varDelta t}}{{\varDelta {z^2}}}{k^2} - {\gamma^1}\frac{{\varDelta t}}{{\varDelta z}}{k^2}} \right] $$

(50)

$$ GP(1) = \left[ { - 2{\beta^1}\frac{{\varDelta t}}{{\varDelta {z^2}}}{k^1} + {\gamma^1}\frac{{\varDelta t}}{{\varDelta z}}{k^1}} \right] $$

(51)

and

$$ EP\left( {NCELL} \right) = \left[ {1 + {\beta^{{NCELL}}}\frac{{\varDelta t}}{{\varDelta {z^2}}}{k^{{NNODE - 1}}} + 2{\beta^{{NCELL}}}\frac{{\varDelta t}}{{\varDelta {z^2}}}{k^{{NNODE}}} - {\gamma^{{NCELL}}}\frac{{\varDelta t}}{{\varDelta z}}{k^{{NNODE}}} + {\gamma^{{NCELL}}}\frac{{\varDelta t}}{{\varDelta z}}{k^{{NNODE}}}} \right] $$

(52)

$$ GP\,\left( {NCELL} \right) = \left[ { - {\beta^{{NCELL}}}\frac{{\varDelta t}}{{\varDelta {z^2}}}{k^{\mathrm{NNODE}}} + {\gamma^{\mathrm{NCELL}}}\frac{{\varDelta t}}{{\varDelta z}}{k^{\mathrm{NNODE}}}} \right] $$

(53)

The same approach was used to numerically solve the transport equation (Eq. 10):

$$ \int\limits_{{k - 1/2}}^{{k + 1/2}} {\alpha \frac{{\partial C}}{{\partial t}}} {\mathrm{d}}z = \int\limits_{{k - 1/2}}^{{k + 1/2}} {\frac{\partial }{{\partial z}}\left[ {{D_{\mathrm{eg}}} + {\alpha_L}\left| v \right|} \right]} \frac{{\partial C}}{{\partial z}}{\mathrm{d}}z + \int\limits_{{k - 1/2}}^{{k + 1/2}} { - \frac{\partial }{{\partial z}}(vC} ){\mathrm{d}}z $$

(54)

Using an implicit time scheme, one obtains:

$$ C_{{k - 1}}^{{t + 1}}G(k) + C_k^{{t + 1}}E(k) + C_{{k + 1}}^{{t + 1}}F(k) = C_k^t $$

(55)

where

$$ G(k) = \left[ {\frac{{\left( {D_{\mathrm{eg}}^{{k - {{1} \left/ {2} \right.}}} + \alpha_L^k{v^{{k - {{1} \left/ {2} \right.}}}}} \right)\varDelta t}}{{{\alpha^k}\varDelta {z^2}}} - \frac{{{v^{{k - {{1} \left/ {2} \right.}}}}\varDelta t}}{{{\alpha^k}\varDelta z}}} \right] $$

(56)

$$ E(k) = \left[ {1 + \frac{{\left( {D_{\mathrm{eg}}^{{k - {{1} \left/ {2} \right.}}} + \alpha_L^k{v^{{k - {{1} \left/ {2} \right.}}}}} \right)\varDelta t}}{{{\alpha^k}\varDelta {z^2}}} + \frac{{\left( {D_{\mathrm{eg}}^{{k + {{1} \left/ {2} \right.}}} + \alpha_L^k{v^{{k + {{1} \left/ {2} \right.}}}}} \right)\varDelta t}}{{{\alpha^k}\varDelta {z^2}}} + \frac{{{u^{{k + {{1} \left/ {2} \right.}}}}\varDelta t}}{{{\alpha^k}\varDelta z}}} \right] $$

(57)

$$ F(k) = \left[ { - \frac{{\left( {D_{\mathrm{eg}}^{{k + {{1} \left/ {2} \right.}}} + \alpha_L^k{v^{{k + {{1} \left/ {2} \right.}}}}} \right)\varDelta t}}{{{\alpha^k}\varDelta {z^2}}}} \right] $$

(58)

Equation (56) can be written in matrix form as follows

Where

$$ G(1) = \left[ { - 2\frac{{\left( {D_{\mathrm{eg}}^1 + \alpha_L^1{\nu^1}} \right)\varDelta t}}{{{\alpha^1}\varDelta {z^2}}} - \frac{{{\nu^1}\varDelta t}}{{{\alpha^1}\varDelta {z^2}}}} \right] $$

(59)

$$ E(1) = \left[ {1 - 2\frac{{\left( {D_{\mathrm{eg}}^1 + \alpha_L^1{\nu^1}} \right)\varDelta t}}{{{\alpha^1}\varDelta {z^2}}} - \frac{{\left( {D_{\mathrm{eg}}^{2} + \alpha_L^1{\nu^2}} \right)\varDelta t}}{{{\alpha^1}\varDelta {z^2}}} + \frac{{{\nu^2}\varDelta t}}{{{\alpha^1}\varDelta {z^2}}}} \right] $$

(60)

$$ F(1) = \left[ {\frac{{\left( {D_{\mathrm{eg}}^{2} + \alpha_L^1{\nu^2}} \right)\varDelta t}}{{{\alpha^1}\varDelta {z^2}}}} \right] $$

(61)

and

$$ G\left( {\mathrm{NCELL}} \right) = \left[ {\frac{{\left( {D_{\mathrm{eg}}^{\mathrm{NNODE}} + \alpha_L^{\mathrm{NCELL}}{\nu^{\mathrm{NNODE}}}} \right)\varDelta t}}{{{\alpha^{\mathrm{NCELL}}}\varDelta {z^2}}} - \frac{{{\nu^{\mathrm{NNODE}}}\varDelta t}}{{{\alpha^{\mathrm{NCELL}}}\varDelta z}}} \right] $$

(62)

$$ E\left( {\mathrm{NCELL}} \right) = \left[ {1 - \frac{{\left( {D_{\mathrm{eg}}^{\mathrm{NNODE}} + \alpha_L^{\mathrm{NCELL}}{\nu^{\mathrm{NNODE}}}} \right)\varDelta t}}{{{\alpha^{\mathrm{NCELL}}}\varDelta {z^2}}} - 2\frac{{\left( {D_{\mathrm{eg}}^{{{\mathrm{NNODE}} + 1}} + \alpha_L^{\mathrm{NCELL}}{\nu^{{{\mathrm{NNODE}} + 1}}}} \right)\varDelta t}}{{{\alpha^{\mathrm{NCELL}}}\varDelta {z^2}}} + \frac{{{\nu^{{{\mathrm{NNODE}} + 1}}}\varDelta t}}{{{\alpha^{\mathrm{NCELL}}}\varDelta z}}} \right] $$

(63)