Published in Vadose Zone Journal 3:848-857 (2004)
© 2004 Soil Science Society of America
677 S. Segoe Rd., Madison, WI 53711 USA
SPECIAL SECTION: RESEARCH ADVANCES IN VADOSE ZONE HYDROLOGY THROUGH SIMULATIONS WITH THE TOUGH CODES
Coupled Vadose Zone and Atmospheric Surface-Layer Transport of Carbon Dioxide from Geologic Carbon Sequestration Sites
Curtis M. Oldenburga,* and
André J. A. Ungerb
a Earth Sciences Division, Lawrence Berkeley National Laboratory, Berkeley, CA 94720
b Earth Sciences Department, University of Waterloo, Waterloo, ON, Canada N2L 3G1
* Corresponding author (cmoldenburg{at}lbl.gov)
Received 9 April 2004.
 |
ABSTRACT
|
|---|
Geologic CO2 sequestration is being considered as a way to offset fossil fuel–related CO2 emissions to reduce the rate of increase of atmospheric CO2 concentrations. The accumulation of vast quantities of injected CO2 in geologic sequestration sites may entail health and environmental risks from potential leakage and seepage of CO2 into the near-surface environment. We are developing and applying a coupled subsurface and atmospheric surface-layer modeling capability built within the framework of the integral finite difference reservoir simulator TOUGH2. The overall purpose of the modeling studies is to predict CO2 concentration distributions under a variety of seepage scenarios and geologic, hydrologic, and atmospheric conditions. These concentration distributions will provide the basis for determining aboveground and near-surface instrumentation needs for CO2 sequestration monitoring and verification, as well as for assessing health, safety, and environmental risks. A key feature of CO2 is its large density (
= 1.8 kg m–3) relative to air ( = 1.2 kg m–3), a property that may allow small leaks to cause concentrations in air above the occupational exposure limit of 4% in low-lying and enclosed areas such as valleys and basements where dilution rates are low. The approach we take to coupled modeling involves development of T2CA, a TOUGH2 module for modeling the multicomponent transport of water, brine, CO2, gas tracer, and air in the subsurface. For the atmospheric surface-layer advection and dispersion, we use a logarithmic vertical velocity profile to specify constant time-averaged ambient winds, and atmospheric dispersion approaches to model mixing due to eddies and turbulence. Initial simulations with the coupled model suggest that atmospheric dispersion quickly dilutes diffuse CO2 seepage fluxes to negligible concentrations, and that rainfall infiltration can cause CO2 to return to the subsurface as a dissolved component in infiltrating rainwater.
|
INTRODUCTION
|
|---|
GEOLOGIC CO2 SEQUESTRATION is being considered as a way of reducing the rate of increase of atmospheric CO2 concentrations due to combustion of fossil fuels for energy production (Reichle et al., 1999). Once injected into deep geologic formations such as depleted oil and gas reservoirs, deep coal seams, and brine formations saturated with groundwater, CO2 will tend to rise upward by buoyant flow even under supercritical conditions unless trapped by low-permeability structures, by dissolution into groundwater, or by mineralogic reactions (Bachu et al., 1994). Despite numerous secondary trapping processes, there is a risk that CO2 will leak from the target storage formation and migrate upwards to where it can seep out of the ground (Oldenburg and Unger, 2003). During leakage and seepage, a fraction of the leaking CO2 may dissolve in groundwater aquifers or in surface waters, thus impacting these natural resources. From the point of view of human and environmental risk associated with exposure to CO2 from leaking geologic carbon sequestration sites, it is advection and dispersion above the ground surface in the biosphere that are most significant, since this is where the key receptors are located. Yet the advection and dispersion processes occurring in the atmospheric surface layer (also referred to simply as the surface layer, and defined approximately as the bottom one-tenth of the atmospheric boundary layer), will be coupled to subsurface processes since (i) the subsurface is the source of the seeping CO2, (ii) ambient air can flow into and out of the subsurface in response to atmospheric pressure changes, and (iii) CO2 is a dense gas that will tend to migrate downwards and hug the ground relative to ambient air. Therefore, simulation models for atmospheric dispersion of CO2 that neglect processes involving the subsurface and the vadose zone in particular may not be appropriate, except in certain limited situations.
A schematic of potential CO2 leakage and seepage from a geologic sequestration site is shown in Fig. 1
, along with associated processes and features. Specifically, Fig. 1 shows a geologic CO2 sequestration site with a permeable fault through which CO2 is unexpectedly leaking upward by buoyancy and pressure-driven flow. The leaking CO2 plume spreads as it decompresses with rise in the subsurface and eventually flows out from the water table into the vadose zone, where it displaces and mixes with existing soil gas. In the vadose zone, leaking CO2 may spread out and upwards until it seeps out of the ground. Above the ground surface in the surface layer, wind and possibly density-driven flow effects will control the flow and dispersion of CO2. A schematic of an eddy-flux tower and CO2 monitoring vault are shown in the figure to suggest potentially important near-surface monitoring strategies for CO2 seepage detection and CO2 sequestration verification.
View larger version (77K):
[in this window]
[in a new window]
|
Fig. 1. Sketch of unexpected leakage and seepage of CO2 from a geologic C sequestration site showing the subsurface and (atmospheric) surface-layer regions, and eddy flux tower and monitoring vault (not to scale).
|
|
Motivated by the need to predict CO2 concentrations in the event that a geologic sequestration site would leak, leading to significant upward migration through the saturated and vadose zones and eventual CO2 seepage at the ground surface, we have developed a coupled subsurface–surface-layer simulation capability called T2CA (TOUGH2 CO2 and Air). T2CA can be used for risk assessment and for designing instrumentation and strategies for geologic carbon sequestration monitoring and verification. This new simulation capability can be used to answer questions about what the expected concentrations will be in the surface layer and shallow subsurface resulting from assumed leakage fluxes. This information can then be used to (i) assess the potential exposure to CO2 for humans and other environmental receptors and (ii) develop specifications and designs of monitoring equipment and strategies for sequestration verification (e.g., Oldenburg and Lewicki, 2003).
In the case of catastrophic failures involving large seepage fluxes (e.g., due to a well blowout), the health risks are obvious and the event could have potentially lethal effects, thus subordinating the sequestration verification issue relative to safety assurance. However, we expect the challenging issues to be health, safety, and environmental risk assessment, as well as monitoring and verification, associated with diffuse or very slow seepage phenomena that are hard to detect. For this reason, our simulation capability is designed for cases of diffuse seepage as opposed to catastrophic failures. The time scale of interest is from 1 mo to 10 yr, making averaging of winds and other environmental variables defensible. The purpose of this paper is to present our approach to modeling subsurface and surface-layer CO2 migration and dispersion of leakage and seepage from geological C sequestration sites, and to show some initial results. This modeling effort is the subject of ongoing testing and verification.
|
Background
|
|---|
Carbon dioxide is a dense gas ( = 1.8 kg m–3) relative to air ( = 1.2 kg m–3), as shown in Fig. 2
, where we have plotted gas density and viscosity for mixtures of CO2 and air calculated from the NIST14 Database (National Institute of Science and Technology, 1992; Magee et al., 1994). Although CO2 is ubiquitous and essential to life as part of the natural C cycle, it is hazardous at high concentrations. The current ambient CO2 concentration in the atmosphere is approximately 375 ppmv (0.0375%); concentrations of 4% can cause immediate danger to humans (NIOSH, 1981). As such, CO2 can be considered a dense hazardous gas, a class of substances that has received considerable attention for leak and spill risk assessment of industrial gases (e.g., Britter and Griffiths, 1982; Hanna and Steinberg, 2001). For example, liquefied propane gas (LPG), liquefied natural gas (LNG), and many others are dense hazardous gases on release to the atmosphere. Motivated by the need to assess risks associated with the mass production and transport of dense gases, a great deal of experimental, analytical, and modeling work has focused on the problem of dense gas dispersion in the surface layer. This work is summarized in a review article by Britter (1989).
View larger version (30K):
[in this window]
[in a new window]
|
Fig. 2. Mixture density and viscosity at 1 bar in the CO2–air system, showing higher density and lower viscosity of gaseous CO2 relative to air.
|
|
The result of many field experiments of dense gas dispersion processes has been the development of correlations involving the most important parameters controlling atmospheric dispersion, including wind speed, density of released gas, and release flux (Britter and McQuaid, 1988). These correlations were developed based on simple scale and dimensional analysis. One of these correlations relates the seepage flux and average wind speed at an elevation of 10 m to delineate regimes of density-dependent and passive dispersion. In Fig. 3
, we plotted this correlation with values appropriate for CO2–air mixtures for four different length scales (L) of the source area along with the typical amount of CO2 emitted and taken up by plants, soil, and roots known as the net ecosystem exchange (e.g., Baldocchi and Wilson, 2001). The length scale L is a characteristic length scale that describes the size (e.g., diameter) of the seepage source region. As shown in Fig. 3, seepage fluxes have to be quite high (note logarithmic scale) for windy situations for the resulting dispersive mixing process to be density-dependent. Note that wind conditions are averages over a period of 10 min.
View larger version (38K):
[in this window]
[in a new window]
|
Fig. 3. Correlation for density-dependent and passive dispersion in the surface layer as a function of seepage flux and wind velocity for four different characteristic source area length scales (L) (see Britter and McQuaid, 1988).
|
|
In prior work (Oldenburg and Unger, 2003), we simulated subsurface migration of leaking CO2 through the unsaturated zone with rainwater infiltration for various leakage rates specified at the water table. Typical seepage fluxes were on the order of 10–5 to 10–6 kg m–2 s–1. As shown in Fig. 3, seepage fluxes of this magnitude lead to passive (i.e., not density-dependent) dispersion for all but the calmest wind conditions. Therefore, we have developed an approach for the case of diffuse emissions that models passive mixing and does not consider dense-gas dispersion effects nor catastrophic CO2 emission events. For these diffuse gas seepage scenarios, we are considering 10-m to 1-km length scales, and time scales from 1 mo to 10 yr.
|
COUPLED MODELING APPROACH
|
|---|
To simulate the coupled subsurface–surface-layer advection and dispersion of CO2, we have developed T2CA, an extension of the EOS7R module of TOUGH2 (Oldenburg and Pruess, 1995; Pruess et al., 1999). TOUGH2 is an integral finite difference reservoir simulator developed for handling multicomponent and multiphase porous media flow systems such as geothermal and hydrocarbon reservoirs, and multiphase flow systems such as the vadose zone or saturated systems that contain contaminant plumes with nonaqueous phase liquids (Pruess et al., 1999). T2CA handles five components (H2O, brine, CO2, a gas tracer, air) and heat. Real gas mixture properties are calculated so the full range of high-pressure sequestration site conditions to low-pressure ambient surface-layer conditions can be modeled. We have added atmospheric surface-layer dispersion capabilities to T2CA to create a fully coupled subsurface–surface-layer simulator. The advantages of adding atmospheric dispersion process-modeling capabilities to the reservoir simulation code as opposed to coupling an existing atmospheric dispersion code to the reservoir code are: (i) consistent multiphase and multicomponent treatment in the subsurface and surface layer for convenient mass-conservative transport, (ii) full coupling of multiphase and multicomponent flow and transport between the subsurface and surface-layer regions, (iii) synchronous time-stepping in the two regions, (iv) lack of need to license or purchase externally developed software, and (v) expediency for us due to our long experience with TOUGH2. For other surface-layer modeling objectives, such as modeling dispersion over nonflat terrain, density-dependent flow, or high Reynolds number flows, coupling of existing Navier–Stokes codes to TOUGH2 is likely to be the most expedient approach.
The purpose of this section is to present the methods implemented in T2CA. Because subsurface transport in T2CA is unchanged from the standard approach used in TOUGH2, we focus our discussion on the methods we apply in the surface layer to model atmospheric dispersion. These methods are derived from the atmospheric dispersion modeling literature and transferred into the TOUGH2 reservoir simulation framework in an expedient way. A simple verification problem is presented to show that the methods are implemented correctly. While the discussion below focuses on CO2 transport, all of the gas-phase components are treated identically.
Transport of Dilute CO2 as a Passive Gas
Transport of CO2 as a passive gas implies that it advects and disperses in the atmosphere without influencing the flow field. In order for this assumption to hold, CO2 must be at concentrations sufficiently low so that it does not significantly affect the density or viscosity of the ambient atmosphere. Under this assumption, we discuss below the underpinnings of the use of an ambient wind profile as well as advection and dispersion processes in the lower layers of the atmosphere as developed in the atmospheric transport literature (e.g., Slade, 1968; Pasquill, 1974; Stull, 1988; Arya, 1999).
Logarithmic Velocity Profile
The ambient time-averaged wind profile near the ground surface has been shown theoretically to follow a logarithmic profile. An excellent review of the assumptions and calculations involved in the logarithmic profile, as well as experimentally derived parameters obtained from calibration to field data, is provided in Slade (1968)(p. 73). The logarithmic wind profile is valid over approximately the lower one-tenth of the atmospheric boundary layer, or approximately a few tens of meters above the ground surface. The logarithmic wind profile as shown on Fig. 4
is given as
 | [1] |
where u(z) is the ambient wind velocity as a function of height, u* is the friction velocity (a parameter that governs the shape of the wind profile near the ground surface for various surface types), k is von Karman's constant (k = 0.4), z is the elevation, and z0 is a roughness height such that u(z) = 0 for z 
z0 and is also a function of various surface types (Slade, 1968). The logarithmic wind profile is strictly applicable only to neutral stability conditions, although equations that account for its variation with atmospheric stability can be formulated (e.g., Golder, 1972).
View larger version (18K):
[in this window]
[in a new window]
|
Fig. 4. Schematic of the logarithmic velocity profile used to approximate time-averaged winds in the surface layer.
|
|
Advective–Dispersive Transport
Gradient Transport Theory
The mean turbulent transport of CO2 as a passive gas in the surface layer can be described by the advective–dispersive transport equation with variable eddy diffusivities (Kx, Ky, Kz) (Arya, 1999, p. 137). For the three-dimensional (x, y, z) transport of a component (such as CO2) at concentration c, this equation is
 | [2] |
For convenience in surface-layer transport modeling, the coordinate system can be arranged so that x is aligned in the downwind direction, making v = w = 0, where u is the ambient wind.
Gaussian Plume Model
For the special case of constant eddy diffusivities and the assumption of a uniform wind velocity (u) with no shear (i.e., no velocity gradient), and assuming that advection dominates diffusion in the x direction, solutions to Eq. [2] are given by the well-known Gaussian plume dispersion model, with eddy diffusivities Dxx, Dyy, and Dzz given by
 | [3] |
where 
x, y, z are the standard deviations of concentration distributions at an observation or receptor point, and t is the travel time to the point (e.g., Arya, 1999, p. 132).
The fundamental challenge in Gaussian plume modeling is the estimation of the eddy diffusivities. The empirically derived Pasquill–Gifford (P–G) dispersion curves provide a practical means of determining atmospheric dispersion and are discussed in detail in Slade (1968) and Arya (1999). Essentially, eddies that are smaller than the plume size are assumed to result in dispersion of passive constituents that can be mathematically represented as a diffusion process. The Pasquill–Gifford dispersion curves were developed from experiments conducted over a wide variety of terrain, such as the project Prairie Grass and British diffusion experiments (Pasquill, 1961; Gifford, 1961), and a variety of atmospheric conditions (Class A, extremely unstable; Class B, moderately unstable; Class C, slightly unstable; Class D, neutral; Class E, slightly stable; Class F, moderately stable). The Pasquill–Gifford curves provide values of y and z as a function of a downwind observation or receptor location under a specific atmospheric condition (Classes A–F) from which constant values of Dyy and Dzz can be derived from Eq. [3]. The Pasquill–Gifford dispersivities are valid for dispersion over distances less than approximately 1 km downwind from near-surface sources over moderately rough and flat terrain (Arya, 1999, p. 203).
Despite the agreement with field data and widespread acceptance of large-scale modeling, the Gaussian plume model assumes uniform velocity, which is not valid in the surface layer near the ground surface, an area of particular interest for CO2 leakage and seepage studies. Before presenting the approach that we have used in this study, we present a verification problem that makes use of the simple analytical solutions of the Gaussian plume model to confirm our implementation of eddy diffusivities and velocity specification in T2CA.
Verification
The simple Gaussian plume model is useful for verifying the surface-layer dispersion capabilities we have developed in T2CA. For a three-dimensional system with uniform wind speed (u) of 1 m s–1, Dxx = Dyy = Dzz = 5 m2 s–1, and a point source strength of 0.314 kg s–1, we solved the Gaussian plume analytical solutions given in Eq. 6.38 and 6.42 of Arya (1999) and compared them with the T2CA solution of the same problem. Note that we take advantage of symmetry in the problem and carry out the T2CA simulation in the upper half of one side of the three-dimensional domain, and use a corresponding source strength that is one-fourth that used in the analytical solution. We show in Fig. 5
T2CA results for the three-dimensional CO2 concentration plume as it is advected by a uniform wind in the x direction (u = 1 m s–1) from a 1 by 1m source with strength Q = 0.0785 kg s–1 in a finely discretized region near the origin (x, y, z < 10 m) and disperses equally in the y and z directions. Note that the CO2 concentration shown in Fig. 5 is in units of kilograms CO2 per cubic meter of gas to match the units of the analytical solutions. Furthermore, the analytical solution assumes isotropic dispersion, whereas in T2CA we assume that advection dominates dispersion in the x direction. To make up for this difference, we used a grid with 10-m gridblocks (
x = 10 m) throughout most of the domain (x, y, z > 10 m) to make numerical dispersion in T2CA approximately match the Dyy and Dzz of the analytical solution, where numerical dispersion in the upstream-weighted and implicit T2CA is approximately x/2 x u = 5 m x 1 m s–1 = 5 m2 s–1.
View larger version (75K):
[in this window]
[in a new window]
|
Fig. 5. T2CA results of CO2 concentration (kg CO2 m–3 gas) for the three-dimensional Gaussian plume verification problem.
|
|
In Fig. 6
we show contours of CO2 concentration for the y–z plane extracted from the three-dimensional domains of both the analytical and T2CA results. As shown, the agreement is very good. For x = 100 m, we have extracted the profile in the y direction and plotted CO2 concentration against y for the analytical and numerical T2CA results as shown in Fig. 7
. As shown, the agreement is very good, and the calculated standard deviations match closely (Fig. 7). Using Eq. [3] with Dyy = Dzz = 5 m2 s–1, the theoretical standard deviation at x = 100 m would be 31.6 m, in good agreement with calculated results. The 10% puff-radius approximation (Arya, 1999, p. 132) matches the theoretical result to within 1% (rp = 67.4 m/2.15 = 31.3 m). These results serve to verify the atmospheric dispersion framework built into T2CA.
View larger version (24K):
[in this window]
[in a new window]
|
Fig. 6. Comparison of T2CA results against analytical solutions for the three-dimensional Gaussian plume verification problem for the x–y plane.
|
|
View larger version (28K):
[in this window]
[in a new window]
|
Fig. 7. Comparison between T2CA results and analytical solution of the CO2 concentration profile in the y direction, with calculated standard deviations.
|
|
Variable-K Theory
Although attractive for its simplicity and widely used, the Gaussian plume model is not valid for situations with wind shear (i.e., a nonzero gradient of u with height), as appropriate for winds near the ground surface that will affect CO2 seepage (Arya, 1999, p. 197–199). Instead, theory and data identify the need for variable eddy diffusivities (Kx, Ky, Kz), an approach called variable-K theory. The variable-K theory is recommended for cases with wind shear and nonhomogeneous turbulence such as will be found in the surface layer (Arya, 1999, p. 143). For our surface-layer applications involving CO2 seepage, we used variable-K theory and the assumption that Kz increases linearly with height as
 | [4] |
(Arya, 1999, p. 143). This model assumes neutral stability in the surface layer, allows for a variable wind speed with height, and models the larger dispersion that occurs as the plume moves upward. There is no analogous formulation of Ky valid for short travel distances (<10 km) in variable-K theory (Arya, 1999, p. 151). Because of this shortcoming of the theory and the urgent need to understand potential leakage and seepage CO2 concentrations, here we adopt a two-dimensional configuration for our test problem that models only vertical dispersion and downwind advection by wind with a logarithmic velocity profile. Because CO2 dispersion will occur only in the vertical direction, this represents a conservative model in that actual CO2 concentrations downwind will be lower for emissions from any realistic areal source for which lateral dispersion occurs. The neglect of lateral dispersion is not an inherent limitation of T2CA, which is in fact three-dimensional, and can include lateral dispersion assuming a reasonable parameterization is available.
Summary
We have combined the logarithmic velocity profile and variable-K theory into a preliminary and expedient approach for modeling multicomponent (CO2, gas tracer, and air) transport in a two-dimensional surface layer that is directly coupled with a porous medium subsurface region. In this approach, we calculate eddy diffusivities from the variable-K diffusivity of Eq. [4] to produce an effective atmospheric dispersivity at every gridblock in the surface layer, a convenient approach in the discretized framework of T2CA. Although it is normally negligible, the molecular diffusion coefficient is added to the eddy diffusivity with the largest term controlling the dispersion process. The single effective dispersivity is then used in the advective–dispersive transport equation for each chemical component to model surface-layer transport. The methods implemented in T2CA for surface-layer dispersion are the subject of ongoing verification and testing.
|
IMPLEMENTATION IN TOUGH2
|
|---|
Specification of the Logarithmic Wind Profile
The simulation of atmospheric advection and dispersion by the above methods requires the specification of a logarithmic wind profile within the TOUGH2 framework that will prevail throughout the simulation. This step involves generating a grid with sufficient layers (i.e., parallel to the ground surface) to discretize the wind profile to the desired accuracy. Next, a static gas-phase pressure profile in the z direction is used along with a constant pressure difference between the upstream and downstream boundaries of the surface layer:
 | [5] |
where P1 and P2 are the upstream and downstream pressures, respectively, within a layer. TOUGH2 computes the phase velocity using Darcy's equation:
 | [6] |
where kD is the intrinsic (Darcy) permeability, 
is the porosity, µ is the gas viscosity, is the mass density of the gas phase, g is the gravitational acceleration, and z is height. Setting the porosity of the surface-layer materials to unity, the velocity of the atmospheric air will be proportional to the permeability of the layer and pressure difference, P, for horizontal layers. Given that P is a constant for all layers, the individual permeability variations of the layers will combine to produce the logarithmic wind profile. Note that the thickness of each layer must be constant to ensure a constant air velocity within the layer across the length of the domain. Note further that the permeability is a pseudopermeability with no physical significance; its purpose is simply to create the desired velocity profile. Note further that the velocity in the surface layer does not change during the simulations because the dispersion process is passive. In essence, we have specified a velocity field for the surface layer that persists throughout the T2CA simulation. In the example presented below, the permeability for the top atmospheric layer with the highest (reference) velocity is set to 1 x 10–2 m2 to minimize P in Eq. [5] and the potential for artificial forced flow of atmospheric air from the upstream boundary into the subsurface. In addition, this value also permits smooth convergence of the Newton iteration.
Calculating Atmospheric Dispersion
Within the TOUGH2 framework, transport of CO2 as a passive gas will follow the advective–dispersive transport equations used to calculate the multicomponent transport of species in the gas phase. Ambient atmospheric dispersion of CO2 is implemented by using a spatially dependent effective molecular diffusivity in the surface-layer region. With this approach, the diagonal of the tensor representing diffusion of CO2 is modified to be the sum of the eddy diffusivity and molecular diffusion.
Numerical dispersion in the implicit and upstream-weighted TOUGH2 framework is on the order of one-half the grid spacing multiplied by the velocity. Because of the alignment of the grid with the unidirectional wind, numerical dispersion occurs only in the flow direction (i.e., x direction) in the surface layer. In the quasisteady cases we are considering, advection dominates transport in the flow direction. In the vertical direction, the velocity is zero (w = 0); thus, vertical eddy diffusion is untainted by numerical dispersion. If CO2 front tracking in the surface layer ever arises as a focus of interest, special weighting schemes can be implemented to diminish numerical dispersion in the flow direction (e.g., Oldenburg and Pruess, 2000).
Restriction to Passive Dispersive Transport
In general, CO2 dispersion can occur both as a dense or as a passive gas, depending on CO2 concentration. Although our approach is applicable only to passive gas transport in the surface layer, note in Eq. [6] that the body force term remains. Therefore, if significant density effects ever arise in the surface layer, velocity will be affected and will deviate from the logarithmic velocity profile that should remain unaltered throughout the simulation. If the velocity profile in the surface layer changes, it is an indication that the atmospheric dispersion process is not strictly passive, and the user should proceed carefully to assess whether other methods should be applied to model dense gas dispersion. Full density dependence is assumed in the subsurface (porous medium) regions, where CO2 concentrations can be quite large and density-driven flow correspondingly important.
Summary
Implementing the coupled subsurface–surface-layer CO2 flow and transport model in TOUGH2 involves the assumption of an average logarithmic wind velocity profile and the use of an effective dispersivity formed by summing the eddy diffusion and molecular diffusion coefficients. Our approach is novel in that it implicitly couples the surface layer to the subsurface region. This coupling is important because CO2 seepage may return to the subsurface through gas-phase advection, diffusion, or dissolution in infiltrating water. While our multicomponent transport methods for the subsurface are firmly established and accepted, we present our surface-layer transport and dispersion approach as a preliminary and expedient multicomponent method useful for estimating surface-layer CO2 concentrations resulting from CO2 leakage.
|
PRELIMINARY RESULTS
|
|---|
We present preliminary results to demonstrate the capabilities of T2CA. The properties of an idealized two-dimensional unsaturated zone and atmospheric surface layer are shown in Fig. 8
, with properties given in Table 1. The domain consists of 100 x 1 x 55 gridblocks in the x, y, and z directions. The subsurface consists of 35 layers of gridblocks of dimension z = 1.0 m. The atmospheric surface layer consists of 20 layers of gridblocks of dimension z = 0.5 m. The horizontal discretization is x = 10 m and y = 1.0 m and is uniform throughout the domain. The bottom boundary is held at constant pressure, while the top boundary is closed. The side boundaries in the subsurface are closed, while the side boundaries in the surface layer are held at constant pressure to generate the logarithmic velocity profile.
In the model system, CO2 is being injected at the water table to model the arrival of leaking CO2 from a deep geologic sequestration site. The CO2 migrates upwards through the unsaturated zone and seeps out of the subsurface into the surface layer. We inject pure water at a constant rate of 10 cm yr–1 at the ground surface to model rainfall infiltration. This rainfall infiltration is capable of transporting dissolved CO2 from the surface layer back into the subsurface, as will be shown below. The subsurface part of this system is a Cartesian version of the radial system we have studied earlier (Oldenburg and Unger, 2003). We use the same leakage rate of 0.1% yr–1 for an assumed 4 x 109 kg CO2 sequestration site, giving rise to a leakage rate of 4 x 106 kg yr–1. If we assume this leakage occurs over 104 m2, the seepage flux is approximately 1.3 x 10–5 kg m–2 s–1. Here we assume a two-dimensional system with no lateral dispersion (Dyy = Ky = 0), and we assume a closed top boundary, both of which cause CO2 concentrations to be larger than in a three-dimensional system with a thicker surface layer. The neglect of lateral dispersion and 10-m surface-layer height are consequences of the choice of test problem and not inherent limitations of T2CA, which is in fact three-dimensional, with no limits on domain height.
The surface-layer part of the system has porosity equal to unity and a logarithmic velocity profile for neutral stability conditions that we specify by using variable permeabilities in the layers above the ground surface, as described above in Specification of the Logarithmic Wind Profile. We defined a reference velocity at an elevation of 10 m above the ground to be 1 and 5 m s–1 to test two different wind conditions. The simulation was run for 6 mo, allowing time for the CO2 to migrate upward through the unsaturated zone and seep out of the ground, where it is advected and dispersed by wind in the atmospheric surface layer. The simulation is isothermal at 15°C.
Results of CO2 mass fraction in the gas phase are shown in Fig. 9a and 9b for wind velocities of 1 and 5 m s–1. As shown in the figures, concentrations of CO2 are quite high in the unsaturated zone because the CO2 sweeps through the pores and displaces existing soil gas with little chance for attenuation (Oldenburg and Unger, 2003). A sharp gradient in concentration is maintained at the ground surface because of the large amount of dilution afforded by the wind, which advects air into the seeping CO2 and carries it downwind. Note that we assumed zero background CO2 concentration in the system to examine the CO2 added by the leakage and seepage processes. As shown in Fig. 9a and 9b, the CO2 concentrations rise strongly in the subsurface, but the CO2 concentrations in the surface layer due to this seepage flux and wind condition are practically negligible. Indeed, Fig. 9a shows that the concentrations increase by approximately 0.0001 by mass fraction (
66 ppmv) just above the source area, and by far less several meters above and downwind from it. Such concentration increases would be easily detectable relative to a background CO2 concentration of 375 ppmv (5.7 x 10–4 mass fraction), but would not be a health hazard (NIOSH, 1981). Dispersion is higher in the 5 m s–1 case than in the 1 m s–1 case because Kz increases with friction velocity and because of the wind dilution effect. The concentrations in the surface layer are essentially steady by t = 6 mo, whereas the concentrations in the subsurface associated with the downward infiltration of rainwater containing dissolved CO2 are still evolving.
View larger version (80K):
[in this window]
[in a new window]
|
Fig. 9. Gas phase mass fraction of CO2 and gas velocity in the coupled subsurface–surface-layer model domain 6 mo after CO2 seepage begins for reference velocity of (a) u = 1 m s–1, and (b) u = 5 m s–1. (c) Liquid saturation and (d) mass fraction of CO2 in the liquid with water velocity for infiltration of 10 cm yr–1 (largest water velocity vector 2 x 10–8 m s–1).
|
|
Note also in Fig. 9a and 9b the apparent subsurface dispersion of CO2 to the right (downwind) of the main subsurface plume. Some of this CO2 is re-entering the subsurface as a dissolved component in infiltrating rainwater. The infiltration source is in the first row of subsurface gridblocks, which obtain CO2 from the surface-layer plume by gas-phase diffusion. Although infiltration in the model is pure water, natural infiltrating rainwater does have significant capacity to dissolve additional CO2 relative to its CO2 content when in equilibrium with ambient atmosphere. Specifically, water in equilibrium with air with 375 ppmv CO2 would contain approximately 0.6 mg CO2 L–1, whereas the solubility of CO2 in water at ground-surface conditions is approximately 1500 mg L–1. Thus, rainwater can dissolve additional CO2 from high-concentration leakage or seepage plumes and transport CO2 downward as a dissolved component. The process of downward reflux of CO2 by water infiltration points out the need for coupled modeling approaches that include interactions between the surface layer and subsurface that may be significant in some situations.
Figures 9c and 9d show liquid saturation and mass fraction CO2 in the liquid phase, respectively. These results show the multiphase and multicomponent aspects of the model inherent to the TOUGH2 framework. Note the downward infiltration that occurs and the attenuating effect of CO2 solubility in water infiltrating into the vadose zone.
Figure 10
shows the CO2 gas-phase mass fractions at a receptor located at the ground surface at x = 645 m (100 m downstream from the source) for three different reference wind speeds where a CO2 mass fraction of 10–4 is approximately 66 ppmv CO2. Once again, these results demonstrate that dispersion increases with wind speed, resulting in lower receptor concentrations of CO2. This is a conservative estimate because actual areal sources with lateral dispersion would result in even lower CO2 concentrations for the same seepage flux. Although the results presented here are two-dimensional, T2CA is a fully three-dimensional model, although wind is required to be unidirectional in the x direction.
View larger version (18K):
[in this window]
[in a new window]
|
Fig. 10. Mass fraction of CO2 in the gas phase at a receptor located on the ground approximately 100 m downstream from the source (x = 645 m) for three reference wind velocities.
|
|
|
CONCLUSIONS
|
|---|
We developed a simulation capability for coupled vadose zone and atmospheric surface-layer advection and dispersion of CO2 that may potentially seep from the ground after leaking from geologic carbon sequestration sites. The purpose of such simulations is to provide input to health, safety, and environmental risk assessments, as well as to make specifications for instrumentation needs, and to design monitoring strategies that can be used to verify carbon sequestration and ensure minimal health and environmental risk. The approach we have taken for the dense gas CO2 is to focus on the difficult-to-detect cases of diffuse gas seepage where fluxes are small and surface-layer concentrations are low. In these scenarios, dispersion in the atmospheric surface layer is passive, and the steady logarithmic velocity profile can be used to approximate time-averaged winds under conditions of neutral stability. Variable-K theory is used to estimate atmospheric dispersion in T2CA.
Preliminary application of the method to a two-dimensional CO2 leakage and seepage scenario shows that although high concentrations of CO2 can develop in the subsurface, dispersion strongly attenuates the seepage plume in the surface layer. Our preliminary simulation shows that while such seepage would be readily detectable by conventional instrumentation, which can detect in the parts per million by volume range, the additional CO2 would not constitute a significant health or environmental hazard for the conditions studied. As testimony to the need for coupled models, we observed that infiltration is capable of bringing CO2 back into the subsurface through dissolution into rainwater infiltrating into the subsurface.
|
ACKNOWLEDGMENTS
|
|---|
We thank Philip Maul (Quintessa Ltd.) and an anonymous reviewer for critical reviews. Two anonymous reviewers of an earlier version helped us to refine our methods. We also thank internal LBNL reviewers Bill Riley, Norm Miller, and Marcelo Lippmann for constructive comments. We have also benefited from stimulating discussions with Bill Riley, Sally Benson, Tom McKone, Karsten Pruess, and Robert Hepple (all at LBNL), and Robert Harley (UC Berkeley). This work was supported in part by the Office of Science, U.S. Department of Energy under contract DE-AC03-76SF00098, and by Cooperative Research and Development Agreement (CRADA) between BP Corporation North America, as part of the CO2 Capture Project (CCP) of the Joint Industry Program (JIP), and the U.S. Department of Energy (DOE) through the National Energy Technologies Laboratory (NETL).
|
REFERENCES
|
|---|
- Arya, S.P. 1999. Air pollution meteorology and dispersion. Oxford University Press, New York.
- Bachu, S., W.D. Gunter, and E.H. Perkins. 1994. Aquifer disposal of CO2–Hydrodynamic and mineral trapping. Energy Convers. Manage. 35:269–279.
- Baldocchi, D.D., and K.B. Wilson. 2001. Modeling CO2 and water vapor exchange of a temperate broadleaved forest across hourly to decadal time scales. Ecol. Modell. 142:155–184.
- Britter, R.E. 1989. Atmospheric dispersion of dense gases. Ann. Rev. Fluid Mech. 21:317–344.[Web of Science]
- Britter, R.E., and R.F. Griffiths (ed.) 1982. Dense gas dispersion. Chem. Eng. Monogr. 16. Elsevier, New York.
- Britter, R.E., and J. McQuaid. 1988. Workbook on the dispersion of dense gases. HSE Contract Research Rep. 17/1988. Health Saf. Exec. Rep., Sheffield, UK.
- Gifford, F.A., Jr. 1961. Use of routine meteorological observations for estimating atmospheric dispersions. Nucl. Saf. 2:47–51.
- Golder, D. 1972. Relations among stability parameters in the surface layer. Boundary-Layer Meteorol. 3:47–58.
- Hanna, S.R., and K.W. Steinberg. 2001. Overview of petroleum environmental research forum (PERF) dense gas dispersion modeling project. Atmos. Environ. 35:2223–2229.
- Magee, J.W., J.A. Howley, and J.F. Ely. 1994. A predictive model for the thermophysical properties of carbon dioxide rich mixtures. Res. Rep. RR-136. Gas Processors Assoc., Tulsa OK.
- National Institute of Science and Technology. 1992. NIST database 14 mixture property database. Version 9.08. U.S. Department of Commerce, Washington, DC.
- NIOSH. 1981. Occupational health guidelines for chemical hazards. NIOSH Publ. 81-123. Available at www.cdc.gov/niosh/pdfs/0103.pdf (verified 30 June 2004). U.S. Gov. Print. Office, Washington, DC.
- Oldenburg, C.M., and J.L. Lewicki. 2003. Near-surface monitoring strategies for geologic carbon dioxide storage verification. Rep. LBNL-54089. Lawrence Berkeley Natl. Lab., Berkeley, CA.
- Oldenburg, C.M., and K. Pruess. 1995. EOS7R: Radionuclide transport for TOUGH2. Rep. LBNL-34868. Lawrence Berkeley Natl. Lab., Berkeley, CA.
- Oldenburg, C.M., and K. Pruess. 2000. Simulation of propagating fronts in geothermal reservoirs with the implicit Leonard total variation diminishing scheme. Geothermics 29:1–25.[Web of Science]
- Oldenburg, C.M., and A.J.A. Unger. 2003. On leakage and seepage from geologic carbon sequestration sites: Unsaturated zone attenuation. Available at www.vadosezonejournal.org. Vadose Zone J. 2:287–296.[Abstract/Free Full Text]
- Pasquill, F. 1961. The estimation of the dispersion of windborne material. Meterol. Mag. 90:33–49.
- Pasquill, F. 1974. Atmospheric diffusion. 2nd ed. John Wiley and Sons, Chichester, UK.
- Pruess, K., C. Oldenburg, and G. Moridis. 1999. TOUGH2 user's guide. Version 2.0. Rep. LBNL-43134. Lawrence Berkeley Natl. Lab., Berkeley, CA.
- Reichle, D., et al. 1999. Carbon sequestration research and development. DOE/SC/FE-1. USDOE, Washington, DC.
- Slade, D.H. (ed.) 1968. Meteorology and atomic energy. U.S. Atomic Energy Commission, Division of Technical Information, Washington, DC.
- Stull, R.B. 1988. An introduction to boundary layer meteorology. Kluwer Academic, Dordrecht, The Netherlands.
- van Genuchten, M.Th. 1980. A closed-form equation for predicting the hydraulic conductivity of unsaturated soils. Soil Sci. Soc. Am. J. 44:892–898.[Abstract/Free Full Text]
This article has been cited by other articles:
|
|
|
|
 
G. Schnaar and D. C. Digiulio
Computational Modeling of the Geologic Sequestration of Carbon Dioxide
Vadose Zone J.,
May 21, 2009;
8(2):
389 - 403.
[Abstract]
[Full Text]
[PDF]
|
|
|
TOP
ABSTRACT
INTRODUCTION
