Dense gas dispersion modeling of CO2 released from carbon capture and storage infrastructure into a complex environment

International Journal of Greenhouse Gas Control 17 (2013) 127–139

Contents lists available at SciVerse ScienceDirect

International Journal of Greenhouse Gas Control journal homepage: www.elsevier.com/locate/ijggc

Dense gas dispersion modeling of CO2 released from carbon capture and storage infrastructure into a complex environment Kun-Jung Hsieh a , Fue-Sang Lien b , Eugene Yee c,∗ a b c

Waterloo CFD Engineering Consulting Inc., Waterloo, ON, N2T 2N7, Canada University of Waterloo, Waterloo, ON, N2L 3G1, Canada Defence Research and Development Canada Suffield, P.O. Box 4000, Medicine Hat, AB, T1A 8K6, Canada

a r t i c l e

i n f o

Article history: Received 22 December 2012 Received in revised form 11 April 2013 Accepted 9 May 2013 Keywords: Atmospheric dispersion modeling Boussinesq approximation Carbon capture and storage Dense gas dispersion Risk assessment

a b s t r a c t Two scenarios for atmospheric dispersion relevant for consequence assessment associated with the loss of containment from carbon capture and storage related infrastructure were investigated using a physicsbased mathematical model: namely, the leakage of carbon dioxide (CO2 ), which is a heavier-than-air (or, dense) gas, from storage tanks and transportation pipelines. Simulations of these two scenarios (viz., a storage tank release in the vicinity of a cubical obstacle and a pipeline rupture in a complex topography involving two axisymmetric hills) were performed using computational fluid dynamics, in which the density variations of the fluid (containing the dense gas) were simplified using the Boussinesq approximation. It is shown that the presence of an obstacle and/or complex terrain has a significant influence on the dispersion of the dense gas. Owing to the ‘slumping’ of the dense gas under the action of gravity, regions well upwind of the source of the gas release can also lie within the hazard zone. The research reported herein provides an improved model for analyzing hazards associated with the dispersion of dense gas clouds and their interaction with local building wakes and/or topographic (terrain) features and contributes to providing a sophisticated method for the assessment of safety and security related to the transportation and geological storage of CO2 . © 2013 Elsevier Ltd. All rights reserved.

1. Introduction Carbon capture and storage (CCS) is a process that involves the capture and storage of carbon dioxide (CO2 ) emissions in deep geological formations. This technology is now widely accepted as a viable method for reducing significantly greenhouse gas emissions into the atmosphere and is being pursued and implemented globally. Carbon dioxide emitted from sources such as power plants and other industrial processes is compressed and transported by pipelines and/or tankers to storage sites, and then injected deep underground into various types of geological structures and formations (e.g., depleted oil and gas reservoirs, deep saline aquifers, unmineable coal seams) for long-term storage.

Abbreviations: CCS, Carbon Capture and Storage; CDS, Central Differencing Scheme; CFD, Computational Fluid Dynamics; ERCOFTAC, European Research Community on Flow, Turbulence and Combustion; IAHR, International Association for Hydro-Environment Engineering and Research; LES, Large Eddy Simulation; NIOSH, National Institute for Occupational Safety and Health; RANS, Reynolds-Averaged Navier–Stokes; SIMPLE, Semi-Implicit Method for Pressure-Linked Equations; SGS, Subgrid Stress; UMIST, Upstream Monotonic Interpolation for Scalar Transport. ∗ Corresponding author. Tel.: +1 403 544 4605; fax: +1 403 544 3388. E-mail addresses: [email protected] (K.-J. Hsieh), [email protected] (F.-S. Lien), [email protected], [email protected] (E. Yee). 1750-5836/$ – see front matter © 2013 Elsevier Ltd. All rights reserved. http://dx.doi.org/10.1016/j.ijggc.2013.05.003

A significant amount of research has been conducted on the geological storage of CO2 . Hosa et al. (2011) and Michael et al. (2010) reviewed various existing and completed CCS projects worldwide that involved the injection of CO2 into saline formations (e.g., Sleipner project in Norway). The operational and reservoir characteristics of injection sites, such as the rate and depth of injection, and the permeability and porosity of the aquifer, were also summarized in these studies. These projects are useful for demonstrating the feasibility of safe CO2 storage, for verifying the monitoring techniques required for CO2 leakage detection, and for providing the detailed experimental data needed for the calibration and validation of mathematical models for prediction of the various complex physico-chemical phenomena associated with CO2 sequestration. Schnaar and Digiulio (2009) provided a comprehensive review of a number of numerical studies of CCS, including characterization and monitoring of the sequestration site, and assessment of the potential leakage of CO2 through caprock, faults, and abandoned wells to the vadose zone and the overlying atmosphere. Class et al. (2009) and Pruess et al. (2004) compared the predictive performance of various reservoir simulators (e.g., TOUGH2, ECLIPSE, etc.) on several benchmark (albeit hypothetical) problems related to the geological storage of CO2 . The benchmark problems considered by Pruess et al. (2004) were mainly focused

128

K.-J. Hsieh et al. / International Journal of Greenhouse Gas Control 17 (2013) 127–139

Nomenclature c¯ k

mass fraction (concentration) of the kth scalar in the mixture Cε1 , Cε2 , Cε3 closure coefficients in the ε-equation closure coefficient in the specification of the eddy C viscosity t D molecular diffusivity of the scalar acceleration due to gravity g G universal gas constant Gk buoyancy production term H reference height k turbulence kinetic energy Mk molecular weight of species k gauge pressure p¯ P absolute pressure Pk rate of production of turbulence kinetic energy flux Richardson number Rf Rk specific gas constant for species k time t T temperature u¯ i ensemble-mean velocity in the xi -direction u friction velocity turbulent scalar fluxes uj ck ui uj Uref x, y, z y+

Reynolds stresses reference velocity Cartesian coordinates dimensionless wall-normal distance

Greek letters ˛ volumetric expansion coefficient Kronecker delta function ıij ε rate of dissipation of turbulence kinetic energy dynamic molecular viscosity of the fluid   kinematic molecular viscosity of the fluid kinematic eddy (or, turbulent) viscosity t  density of the fluid c turbulent Schmidt number for diffusion of the scalar k turbulent Schmidt number for diffusion of k turbulent Schmidt number for diffusion of ε ε Subscripts i, j, k Cartesian tensor indices 0 reference point

on one-dimensional (1D) homogeneous media, whereas problems involving three-dimensional (3D) heterogeneous media were studied by Class et al. (2009). These included the leakage of CO2 through an abandoned well, the utilization of CO2 injection for enhanced methane recovery, and the injection of CO2 into a heterogeneous (geological) formation. In broad terms, these investigations found that there was generally a good agreement between the predictions provided by the different simulators, and that the discrepancies between these predictions were largely due to different assumptions and simplifications implicit in each simulator. Although a number of investigations have been conducted on the subsurface modeling of CCS, only a few studies have considered the leakage of CO2 from the subsurface to the vadose zone, followed by the release of the CO2 to the atmosphere. Pruess (2008) used TOUGH2 (Pruess et al., 1999) to model the accumulation of CO2 in a secondary storage reservoir (at a depth of about 100 m below the ground surface) resulting from migration of the gas upwards through a fault that connects the reservoir to a deep geological

formation where the CO2 was sequestered. This is followed by the subsequent discharge of the CO2 from the storage reservoir into the atmosphere through still another fault that provides a path from the reservoir to the ground surface above. Ogretim et al. (2009) and Oldenburg and Unger (2004) studied CO2 leakage from the water table through the vadose zone and subsequent dispersion into the atmosphere, whereas de Lary et al. (2012) examined the impacts on human health of CO2 leakage from the vadose zone into a building followed by the subsequent indoor dispersion of the material. Physical phenomena related to the slow seepage of CO2 from the vadose zone to the overlying atmosphere or into a building were investigated in these studies, as opposed to a catastrophic release of CO2 (e.g., arising from a well blowout or a pipeline rupture). Ogretim et al. (2009) examined the effects of crosswind and topographic interaction on the migration of CO2 in the vadose zone under a hill-like topography. The wind velocity over the hill surface was obtained using potential flow theory, and this information was used subsequently to calculate the pressure distribution along the hill surface based on the Bernoulli equation. This pressure distribution was then imposed as the boundary condition at the top of a two-dimensional (2D) subsurface simulation domain, for which the subsurface modeling was conducted using TOUGH2 with the EOS7CA module (the latter of which provided a detailed specification of the water–air–CO2 equations of state in the vadose zone). Along a similar vein, Oldenburg and Unger (2004) developed a coupled subsurface and atmospheric surface layer simulation capability within the framework of TOUGH2. In their coupled modeling approach, the mean transport and turbulent diffusion of CO2 in the atmospheric surface layer over level terrain was treated simply as a passive gas and was modeled using a simple advectiondiffusion equation. The gradient diffusion hypothesis was used as the simplest closure model for the turbulent scalar fluxes in the advection-diffusion equation, and a logarithmic mean wind profile applicable for a neutrally-stratified atmosphere was imposed. Finally, de Lary et al. (2012) simulated the leakage and subsequent migration of CO2 from a single point at the bottom of the vadose zone into a building using TOUGH2 with the EOS7CA module, and subsequently predicted the CO2 concentration in the building using a simple analytical model in order to assess the health impacts on humans due to CO2 leakage from underground. They concluded that the indoor CO2 concentrations only reached hazardous levels under very specific conditions such as a high leakage flow rate from the source and/or a low building ventilation rate. An interesting case study of a real-world incident involving the catastrophic release of CO2 has been conducted recently by Hedlund (2012). In this study, he used numerical simulation to model the CO2 outburst at the Menzengraben potash mine, where several thousand tons of CO2 were blown out of a mine shaft. To this purpose, two source terms associated with high- and low-momentum CO2 releases were modeled separately. Their numerical simulation results showed that in the case of the high-momentum release, the leaked CO2 diluted quickly and never reached the ground surface. In contrast, the low-momentum release resulted in a high concentration gas cloud near the ground surface. Mazzoldi et al. (2011, 2013) conducted quantitative risk assessment of scenarios involving CO2 leakages for various CCS projects using computational fluid dynamics (CFD). In particular, these researchers focused on provision of an improved understanding of elements constituting the risk associated with an accidental release arising from the transport of the captured CO2 from the industrial sources to a suitable storage site along high-pressure CO2 pipelines. In both studies, the pipelines were assumed to be placed above ground level, and the downwind distance from the leakage source reached by a harmful concentration level (100,000 ppmv or larger) of CO2 was predicted and used to assess the risk to human health. Mazzoldi et al. (2011) simulated different release

K.-J. Hsieh et al. / International Journal of Greenhouse Gas Control 17 (2013) 127–139 Table 1 Test cases pertinent to the leakage of CO2 from storage and transportation facilities. Test cases

Release scenario

Thorney Island (Davies and Singh, 1985)

Storage tank leakage near a building Pipeline rupture in the vicinity of hilly terrain

Axisymmetric hill (Simpson et al., 2002)

scenarios for CO2 using both CFD (based on the Reynolds-averaged Navier–Stokes approach) and Gaussian plume models, and evaluated the predictive performance of these two approaches against two field experiments (Mazzoldi et al., 2008). They concluded (perhaps not surprisingly) that CFD models offered improved risk assessments for hazards associated with the dispersion of CO2 clouds compared to the simpler Gaussian plume models. It is important to understand the risks associated with the deployment of CCS, and a principal component of this understanding involves the prediction and quantification of the effects of CO2 leakage from the subsurface where it is stored to the atmosphere where it is dispersed by natural turbulent mixing processes on a variety of scales. The objective of this paper is to investigate the atmospheric dispersion of CO2 arising from its leakage from CCS related infrastructure such as pipelines and CO2 storage tanks. This paper focuses on the physics-based modeling of large and rapid CO2 releases (e.g., from a pipeline rupture or loss of containment in a storage facility) using CFD. These large releases of CO2 will result necessarily in the formation of dense gas clouds (viz., clouds that have greater density than the ambient environment) which have a tendency to disperse at ground level (where people will be subjected to their effects). Indeed, negatively buoyant clouds/plumes of contaminants tend to cover large surface areas with limited vertical extent, leading potentially to more serious environmental consequences because the gas will remain near the ground and populated areas. The dispersion behavior of negatively buoyant releases is complicated even for the simplest conditions (unobstructed and flat terrain), but in reality the complexity of the phenomena is further exacerbated by the interaction of the dense gas cloud with building wakes and/or topographic features. Consequently, the modeling of dense gas dispersion needs to incorporate properly the influences of real terrain and local obstacles on the dispersion. A numerical three-dimensional dense gas dispersion model representing the interaction between the dense gas cloud and obstacles or complex terrain is developed, and applied to the simulation of two scenarios that are similar to (or mimic) the leakages of CO2 from storage and transportation facilities (cf. Table 1). The first test case involves simulation of the Thorney Island Phase II field experiments (Davies and Singh, 1985). This test case can perhaps serve as a prototype for a catastrophic and instantaneous release of a dense gas from a storage tank located near a building, and tests the accuracy of predictions of dense gas dispersion under the influence of a ruptured container and a local building (where the released volume is comparable to the building scale). Alternatively, as pointed out by an anonymous reviewer, this test case can also be a surrogate for a scenario involving the rapid release of a moderate quantity of CO2 gas from a near-surface (perhaps temporary) geological storage site. The second test case involves a simulation of the complex flow over an axisymmetric three-dimensional hill (Simpson et al., 2002), and is one of the test cases used in the 11th European Research Community on Flow Turbulence and Combustion (ERCOFTAC)/International Association for Hydro-Environment Engineering and Research (IAHR) Workshop on Refined Turbulence Modeling. Although there is no dispersion data available for this experiment, the detailed measurements of flow are used to validate our predictions of the mean flow and turbulence fields for the complex topography, and then this information is used subsequently

129

to develop a scenario involving the continuous release of a dense CO2 gas plume (from an accidental or malicious pipeline rupture in a transportation facility) over complex terrain and to investigate numerically the influence of the topographic features on the dense gas dispersion. In both the Thorney Island Phase II field experiments and the axisymmetric three-dimensional hill test case, the Boussinesq approximation for the representation of the buoyancy term in the momentum equation is used to simplify the modeling of the dense gas dispersion. It should be noted that in the case of CO2 storage or transportation facilities (e.g., storage tanks, pipelines), the material is normally stored or transported in a liquefied form after undergoing pressurization. In consequence, it is expected that a high-pressure CO2 release will first experience a significant pressure drop, followed by a substantial temperature drop due to the Joule–Thompson effect. As a result, the CO2 will undergo a phase change from liquid to solid or gas, and potentially result in the formation of dry ice and/or of a sonic jet flow during the process of pressure decompression. Currently, many of the physical processes associated with the rapid release of CO2 have not been explicitly modeled in dense gas dispersion models [see Koornneef et al. (2010) for a discussion of these issues]. Nevertheless, Mazzoldi et al. (2011) argued that most rapid releases of a liquid/vapor mixture of CO2 would transform rapidly into the gas phase due to the frictional heating and mixing provided by the ambient air (following rapid entrainment of air into the initial dense gas cloud). Although source term modeling is an important problem, it is beyond the scope of the current study. Accordingly, we will simply assume [following the arguments by Mazzoldi et al., 2011] that the CO2 liquid/gas mixture has been converted entirely into the gas phase immediately after its release from the source. This is the implicit assumption made in our proposed scenarios described later in Section 3.2. 2. Mathematical formulation and numerical framework The dispersion of a heavier-than-air (dense) gas cloud is predicted by solving the three-dimensional conservation equations for mean (ensemble-averaged) quantities in a turbulent flow field. This involves solving the conservation equations of mass, momentum, and species concentration for an isothermal (constant temperature) fluid flow based on the Reynolds-averaged Navier–Stokes (RANS) approach: ∂ ∂ (u¯ j ) = 0, + ∂t ∂xj

(1)

∂ ∂p¯ ∂ ∂ (u¯ j u¯ i ) = − + (u¯ i ) + ∂t ∂xj ∂xi ∂xj



∂u¯  i ∂xj

 −

∂ (ui uj ) ∂xj

+ ( − 0 )gi ,

∂ ∂ ∂ (¯ck ) + (u¯ j c¯ k ) = ∂t ∂xj ∂xj k = 1, 2, . . . , N,

(2)



∂¯c D k ∂xj

 −

∂ (uj ck ), ∂xj (3)

where u¯ i is the ensemble-mean (or, Reynolds-averaged) velocity in the xi -direction with i = 1, 2, 3 representing the x, y, and z directions, respectively; t is the time; p¯ is the mean (gauge) pressure;  is the dynamic molecular viscosity of the fluid; c¯ k is the mass fraction (concentration) of the kth scalar (or species) in the mixture; D is the molecular diffusivity of the scalar; gi is the gravitational acceleration in the xi -direction;  is the density of the fluid; and, 0 is a reference density (chosen herein as the density of air at sea level with a representative value of 0 = 1.25 kg m−3 ). Summation is implied by repeated indices. Reynolds-averaging of the momentum

130

K.-J. Hsieh et al. / International Journal of Greenhouse Gas Control 17 (2013) 127–139

transport equation (2) and the advection-diffusion equation (3) gives rise, respectively, to the Reynolds stresses which are defined as the tensor ui uj and the turbulent scalar fluxes which are defined as the vector uj ck (k = 1, 2, . . ., N). These equations, along with an equation of state that uses the ideal gas law approximation for the density, namely P

P  = N = , R∗ T T k=1 Rk c¯ k

(4)





 ∂  ∂  ∂   (¯ck − c¯ k,0 ), (5)  = 0 + (P − P ) + (T − T ) + 0 0 ∂P  ∂T  ∂¯ck  0

ı (P − P0 )  ıP   − 0 ≡ = − ˛k (¯ck − c¯ k,0 )≡ − ˛k (¯ck − c¯ k,0 ), 0 0 P0 P0

N

0

k=1

0

with the partial derivatives in Eq. (5) evaluated at the reference state (P0 , T0 , c¯ 1,0 , . . . , c¯ N,0 ) and 0 = (P0 , T0 , c¯ 1,0 , . . . , c¯ N,0 ). For an isothermal atmosphere, the third term on the right-hand side of Eq. (5) vanishes identically. Inserting the ideal gas law of Eq. (4) into Eq. (5), the density perturbations (from the reference state) can be

N

N

k=1

k=1

(6) where ı ≡ ( − 0 ) and ıP ≡ (P − P0 ) = p¯ are used to denote the density and pressure perturbations, respectively; and,



˛k ≡ −

are the main governing equations. Here, P is the absolute pressure, T is the absolute temperature (assumed to be a constant for the isothermal atmosphere considered herein), and Rk ≡ G/Mk is the specific gas constant for species k (G = 8.314472 × 103 J K−1 kg mol−1 is the universal gas constant and Mk is the molecular weight of species k). Also, in Eq. (4), R∗ is the specific gas constant for the mixture consisting of N species. It is noted that the gauge pressure p¯ in Eq. (2) is the pressure deviation from an adiabatic atmosphere in hydrostatic equilibrium with corresponding density0 , with the result that the absolute pressure in Eq. (4) is related to the gauge pressure in Eq. (2) as P = p¯ + P0 , where P0 is some nominal atmospheric pressure level in the adiabatic atmosphere (e.g., at the ground surface). If the density variations are not too large, we can apply the Boussinesq approximation to Eqs. (1)–(3). Indeed, it is expected that the numerical solution of Eqs. (1)–(4) that incorporates the full effect of the density variations on the fluid flow and the mass fraction of the various species is a very challenging problem. In consequence, we consider application of the Boussinesq approximation in order to ensure a more tractable problem. Nevertheless, it is important to stress that the exact limits of validity of the Boussinesq approximation and the precise effects of this simplification on the model predictions are not known quantitatively. To the authors’ best knowledge, this issue has not been addressed by previous investigators and is beyond the scope of the present study. Undoubtedly, a useful future research endeavor would be to consider a detailed comparison of model predictions obtained with and without using the Boussinesq approximation in order to quantify the errors in the predictions arising from use of this approximation. The main step in developing the Boussinesq approximation is to treat the density as a constant (which can be taken as the reference density0 ) in the unsteady, convective and diffusive terms of Eqs. (1)–(3), while retaining the density variations only in the buoyancy (gravity) term ( − 0 )gi in the mean momentum transport equation. This approach was used by Perdikaris and Mayinger (1993) and Scargiali et al. (2005) to simulate dense gas dispersion over a topographically complex terrain. However, the approximation for the buoyancy term given in these two papers is not correctly stated. In consequence, we provide a simple derivation of the Boussinesq approximation of ( − 0 )gi in the context of the dispersion of a heavier-than-air gas released into the atmosphere. Toward this purpose, consider a Taylor series expansion of the density  = (P, T, c¯ 1 , . . . , c¯ N ) to first order about a reference state (denoted using a subscript 0) given by



expressed as follows:

1 ∂  = 0 ∂¯ck 

N

Rk

R c¯ i=1 i i,0

0

,

(7)

is the volumetric expansion coefficient arising from perturbations in the concentration of species k. To simplify the approximation in Eq. (6), we consider order-ofmagnitude estimates for ıP/P0 in comparison to ı/0 . To proceed, it is assumed that the pressure gradient term in the mean momentum transport equation (2) is of the same order of magnitude as the buoyancy term and, with the buoyancy force acting in the y (i = 2) direction, this implies within the Boussinesq approximation that

       1 ∂p¯   1 ∂ıP    − 0   ı   ≡ ≈  g . g ≡   0 ∂y   0 ∂y   0  0

(8)

An estimate of the magnitude of the term on the left-hand-side of Eq. (8) [obtained by using the ideal gas law of Eq. (4) for the reference state] gives

     1 ∂ıP  R∗ T0 ıP   ≈ ,  0 ∂y  P0 ly

(9)

where ly is a representative or characteristic spatial scale of the circulations in the atmospheric surface layer (and typically, ly ∼ 100 m in an adiabatic atmosphere). Insertion of the order-of-magnitude estimate of Eq. (9) into Eq. (8) leads to

  ıP  P0

 

 

 

ly g ı ly g ı ly g0 ı ≈ = = . R∗ T0 0 P0 0 P0 /0 0

(10)

To proceed further, we need an estimate for P0 which can be obtained as follows. The reference pressure P0 at the Earth’s surface is determined for an atmosphere in hydrostatic equilibrium, so P0 ≈ 0 gH0 where H0 is the effective height in the atmosphere where the pressure is equal to zero (approximately). H0 is estimated to be about H0 ∼ 8000 m (Pielke, 2002). Now, using this estimate for P0 in Eq. (10), it follows that

  ıP  P0

 

 

 

ly g0 ı ly g0 ı ly ı ≈ = = , P0 0 0 gH0 0 H0 0

(11)

from which it is seen that (|ıP|/P0 )  (|ı|/0 ) because ly /H0  1. From this observation, it is now evident that the density perturbations in Eq. (6) are well approximated as

  − 0 ı ≡ ≈− ˛k (¯ck − c¯ k,0 ). 0 0 N

(12)

k=1

Finally, for the special case of a two-component mixture consisting of air (a, or species k = 1) and a contaminant species (g, or species k = 2) [CO2 or some other contaminant in the cases considered herein], the density perturbation of Eq. (12) reduces to  − 0 = ˛¯c , 0

(13)

on noting c¯ a = 1 − c¯ g ≡ 1 − c¯ , ˛ ≡ (g − a )/g , and assuming that the concentrations of air and contaminant gas in the reference state are c¯ a,0 = 1 and c¯ g,0 = 1 − c¯ a,0 = 0. Using the density

K.-J. Hsieh et al. / International Journal of Greenhouse Gas Control 17 (2013) 127–139

perturbation given by Eq. (13), the buoyancy term can be approximated as follows:

8

( − 0 )gi = 0 ˛¯c gi .

6

fluxes uj c  , required to close the transport equations for the mean momentum and mean concentration of the contaminant species, are modeled using the Boussinesq eddy-viscosity approximation (which should not be confused with the Boussinesq approximation for buoyant-driven flows) and the simple gradient diffusion hypothesis, respectively, as follows: ui uj

2 = kıij − t 3

uj c  = −



∂u¯ j ∂u¯ i + ∂xi ∂xj

∂ ∂ ∂ (0 ε) + (0 u¯ j ε) = ∂t ∂xj ∂xj



t ∂k 0 k ∂xj





t ∂ε 0 ε ∂xj

+ 0 (Pk + Gk ) − 0 ε, (17)



ε ε2 + Cε1 0 (Pk + Gk )(1 + Cε3 Rf ) − Cε2 0 , k k

(18)

where Pk ≡ −ui uj ∂u¯ i /∂xj is the turbulence kinetic energy production term,  k = 1.0,  ε = 1.3, Cε1 = 1.44, and Cε2 = 1.92 are model (closure) constants. Furthermore, in Eqs. (17) and (18), Gk is the production term relating to buoyancy and is given by Gk = −˛gi

t ∂¯c , c ∂xi

0

5

10

0

5

10

x/H

4 2 0 -2 -4 -6 -10

(16)

where t = C k2 /ε is the kinematic (turbulent) eddy viscosity, C = 0.09 is a model constant,  c = 0.63 is the turbulent Schmidt number for the scalar, k is the turbulence kinetic energy, ε is the dissipation rate of k, and ıij is the Kronecker delta function. The modeled transport equations for k and ε in the standard k–ε model with buoyancy terms assume the following form within the context of the Boussinesq approximation:

-5

6

(15)

t ∂¯c , c ∂xj

∂ ∂ ∂ (0 k) + (0 u¯ j k) = ∂t ∂xj ∂xj

0 -10



,

4 2

z/H

For the simulation of both test case scenarios studied in this paper, Eq. (14) will be used to model the buoyancy term in the mean momentum equation within the framework of the Boussinesq approximation. It is noted that this approximation for the buoyancy term differs from that used in Scargiali et al. (2005) in that the quantity ˛ involves normalization by a rather than g (as used in the current formulation). Finally, for a two-component mixture consisting of air and a contaminant gas, we need only consider a transport equation for the contaminant species concentration c¯ because the concentration of the air in the mixture is determined from (1 − c¯ ). Within the framework of the standard high-Reynolds-number k–ε model, the Reynolds stresses ui uj and the turbulent scalar

y/H

(14)

131

-5

x/H

Fig. 1. A two-dimensional view of the computational mesh for trial number 26 in a vertical x–y plane at z = 0 (top panel) and in a horizontal x–z plane at y = 0 (bottom panel). The mesh consists of 189 × 60 × 83 nodes in the streamwise (x), vertical (y), and spanwise (z) directions, respectively.

scheme (CDS). Advective volume-face fluxes were approximated using a second-order accurate Upstream Monotonic Interpolation for Scalar Transport (UMIST) scheme (Lien and Leschziner, 1994). The transient term was discretized using a fully implicit, secondorder accurate three-time-level method described in Ferziger and Peric (2002). The Semi-Implicit Method for Pressure-Linked Equations (SIMPLE) algorithm (Patankar and Spalding, 1972) was used to ¯ A nonlinear interpolation scheme determine the (gauge) pressure p. (Rhie and Chow, 1983) was used to interpolate the cell face velocities from the adjacent nodal velocities at the cell centers in order to prevent checkerboard oscillations from developing in the pressure field. 3. Results and discussion 3.1. Storage tank release with obstacle (Thorney Island)

(19)

for the case of a two-component mixture (where c¯ is the contaminant species concentration). In Eq. (18), Rf = −Gk /(Pk + Gk ) is the flux Richardson number and Cε3 = 0.8 is an additional model constant (Rodi, 1978) that is applicable both for vertical and horizontal buoyant shear layers. In particular, for a vertical shear layer where the lateral velocity component is normal to the direction of gravity, we note that Rf = 0. This is the case of relevance for our investigation. The system of partial differential equations for atmospheric dispersion modeling was solved numerically using a collocated, finite-volume method and implemented in the urbanSTREAM code (Yee et al., 2007; Lien et al., 2010). Diffusive volume-face fluxes were discretized using a second-order accurate central differencing

Extensive field trials on the dispersion of dense gas clouds at ground level in the atmosphere were performed by the Health and Safety Executive in Great Britain. The trials were conducted on the Thorney Island airfield in southern England from July 1982 to July 1983. In these trials, isothermal gas clouds of fixed volume and of varying density were instantaneously released into the Table 2 Parameters used for the power-law specification of the vertical profile of the streamwise mean wind velocity for Thorney Island trial numbers 26 and 29 (Phase II). Trial number

Wind speed u¯ 0 (m s−1 )

Exponent n

26 29

1.9 5.6

0.07 0.15

132

K.-J. Hsieh et al. / International Journal of Greenhouse Gas Control 17 (2013) 127–139

6

Trial 26: t = 0 s

8

C 100.00 46.42 21.54 10.00 4.64 2.15 1.00 0.46 0.22 0.10

8

C 100.00 46.42 21.54 10.00 4.64 2.15 1.00 0.46 0.22 0.10

8

C 100.00 46.42 21.54 10.00 4.64 2.15 1.00 0.46 0.22 0.10

8

C 100.00 46.42 21.54 10.00 4.64 2.15 1.00 0.46 0.22 0.10

y/H

4

2

-10

-8

-6

-4

6

-2

x/H

0

2

4

6

Trial 26: t = 2.84 s

y/H

4

2

-10

-8

-6

-4

6

-2

x/H

0

2

4

6

Trial 26: t = 14.68 s

y/H

4

2

-10

-8

-6

-4

6

-2

x/H

0

2

4

6

Trial 26: t = 57.32 s

y/H

4

2

-10

-8

-6

-4

-2

x/H

0

2

4

6

Fig. 2. Time evolution of mean flow streamlines and contours of mean concentration (% v/v) of the gas cloud in a vertical plane at z/H = 0 for trial number 26.

atmosphere. The gas concentrations in the cloud, consisting of Freon-12 and air mixtures, were measured using sensors placed at various downwind fetches from the release location. In this paper, we consider the trials conducted in the second phase of the Thorney Island test series (Davies and Singh, 1985). In these experiments, the dense gas cloud dispersion was conducted in the presence of an obstacle. The obstacle used was a sharp-edged solid cube with a dimension of 9 m on each side. In the instantaneous releases, a cylindrical gas tent of 14 m diameter and 13 m height with a total volume of 2000 m3 was filled with a mixture of Freon-12 (dichlorodifluoromethane, CCl2 F2 ) and nitrogen. More specifically, the released gas consisted of a mixture of 68.4% nitrogen and 31.6% Freon-12 (w/w) to give a mixture density that is roughly twice that of air. The top and sides of the cylindrical tent were quickly removed, leaving an upright cylinder of dense gas exposed to the ambient conditions. The mixture immediately began to expand outwards isotropically in all horizontal directions by gravitational slumping. The gravitational slumping and atmospheric dispersion of this dense gas cloud is affected by its interaction with the cubical obstacle, located either downwind or upwind of the release.

In trial number 26, the windward face of the cubical obstacle was located 50 m downwind from the center of the cylindrical gas tent. The measured wind speed at 10 m above the ground was 1.9 m s−1 , and the concentration of the gas cloud was measured on the windward and leeward faces of the obstacle at heights of 6.4 m and 0.4 m above ground level, respectively. In trial number 29, the leeward face of the cubical obstacle was located 27 m upwind from the center of the cylindrical gas tent. The measured wind speed at 10 m above the ground was 5.6 m s−1 , and the concentration was measured on the leeward face of the obstacle at the height of 0.4 m above ground level. Simulations of the flow and dense gas dispersion for trial number 26 were conducted on a computational mesh of 189 × 60 × 83 control volumes (nodes) in a spatial domain (shown in Fig. 1) with an extent of 16.4H × 4.4H × 14.0H in the streamwise, vertical and spanwise directions, respectively, where H = 9 m is the height of the cubical obstacle. The simulations for trial number 29 were performed on a computational mesh of 139 × 60 × 83 nodes in a spatial domain that spanned 15.9H × 4.4H × 14.0H in the streamwise, vertical and spanwise directions, respectively. The mesh used for the computational domain of trial number 29 (not shown) was very

K.-J. Hsieh et al. / International Journal of Greenhouse Gas Control 17 (2013) 127–139

Trial 29: t = 0 s

5

10

C 100.00 46.42 21.54 10.00 4.64 2.15 1.00 0.46 0.22 0.10

10

C 100.00 46.42 21.54 10.00 4.64 2.15 1.00 0.46 0.22 0.10

10

C 100.00 46.42 21.54 10.00 4.64 2.15 1.00 0.46 0.22 0.10

10

C 100.00 46.42 21.54 10.00 4.64 2.15 1.00 0.46 0.22 0.10

y/H

4 3 2 1 0

x/H

5

Trial 29: t = 3.09 s

5

y/H

4 3 2 1 0

x/H

5

Trial 29: t = 8.47 s

5

y/H

4 3 2 1 0

x/H

5

Trial 29: t = 20.52 s

5

y/H

4 3 2 1 0

x/H

5

133

Fig. 3. Time evolution of mean flow streamlines and contours of mean concentration (% v/v) of the gas cloud in a vertical plane at z/H = 0 for trial number 29.

similar to that used in trial number 26, where the grid lines were preferentially concentrated near the solid surfaces in order to better capture the expected sharp gradients of the flow properties here. At the inlet boundary of the spatial domain used for the simulations, the vertical profile of the mean streamwise velocity was n ¯ approximated using a power-law functional form u(y) = u¯ 0 (y/y0 ) , where u¯ 0 is the reference velocity at the reference height y0 = 10 m. Table 2 summarizes the values of u¯ 0 and n used in the simulations for trial numbers 26 and 29. This power-law profile was also used by Qin et al. (2007) and Sklavounos and Rigas (2004) for their simulations of the Thorney Island experiments. At the outlet boundary, zero Neumann boundary conditions were imposed on all flow variables. At the spanwise and top boundaries, symmetry and free-slip boundary conditions were applied to the flow variables, respectively. At all the solid boundaries (ground and obstacle faces), standard wall functions were used for the mean velocities and turbulence quantities, and the zero Neumann boundary condition was applied for the concentration of the gas (implying no flux of material across the solid boundary). In the simulations, computations were conducted firstly to predict the steady-state flow field, including as such the complex

influence of both the cylindrical gas tent and the cubical obstacle on the flow. This disturbed flow field was then used as the initial condition for the subsequent predictions of the dense gas dispersion, which occurred immediately after the collapse of the cylindrical gas tent. The gas mixture inside the cylindrical tent was released instantaneously to the atmosphere at time t = 0 s. Figs. 2 and 3 show the mean flow streamlines and contours of concentration at different times after the dense gas cloud was released and documents the time evolution of the dense gas cloud. The gravity slumping effect on the dispersion can be seen clearly in both figures, particularly at the earlier times where the streamlines within the highly concentrated gas cloud exhibit a strong downward bulk motion. Note that initially the buoyancy generated forces and pressure gradients arising from density differences between the gas cloud and its environment lead to a bulk motion that causes the dense cloud to spread in all directions near the release location (including a lateral spreading, as well as an upwind spreading against the prevailing wind direction). Further downstream from the release, the bulk flow resulting from the gravitational slumping weakens as diffusion, mixing, and entrainment between the gas cloud and the ambient air reduces the negative buoyancy effects of the cloud

134

K.-J. Hsieh et al. / International Journal of Greenhouse Gas Control 17 (2013) 127–139

5

c (% v/v)

4 3 2 1 0

0

50

100

0

50

100

150

200

250

300

150

200

250

300

t (s)

c (% v/v)

3

2

1

0

t (s)

Fig. 4. Time history of gas concentration on the windward face of the obstacle at a height of 6.4 m (top panel) and on the leeward face of the obstacle at a height of 0.4 m (bottom panel) for trial number 26: (––) experimental data from Davies and Singh (1985); (—) simulation.

and strengthens the influence of the externally imposed velocity field on the transport and dispersion of the cloud. It is noted in trial number 26 that the dense gas cloud flows around and over the cubical obstacle located downwind of the release (cf. Fig. 2). In contrast, for trial number 29, the cubical obstacle located upwind of the release is seen to block the further upwind spread of the dense gas cloud arising from the gravitational slumping. Also, in addition to the blocking effect, the gas cloud undergoes increased dilution in the wake of the upstream cubical obstacle (cf. Fig. 3), as well as retention or “hold-up” of the material within the wake of the obstacle followed by a slow detrainment of material from the wake. This retention of material in the wake followed subsequently by the detrainment of the material converts the putative instantaneous release into a time-varying one. Figs. 4 and 5 display a comparison between the measured and predicted gas concentration time histories for trial numbers 26 and 29, respectively. In particular, Fig. 4 shows that the time histories of concentration observed at the windward and leeward faces of 25

c (% v/v)

20 15 10 5 0

0

20

40

t (s)

60

80

100

Fig. 5. Time history of gas concentration on the leeward face of the obstacle at a height of 0.4 m for trial number 29: (––) experimental data from Davies and Singh (1985); (—) simulation.

the cubical obstacle are quite different. Note that the persistence of the gas cloud (as measured by the difference between the arrival and departure times of the cloud) observed at the leeward face of the cubical obstacle is significantly greater than that observed at the windward face. This demonstrates the effect of the obstacle on dispersion, where the recirculation bubble in the wake region of the obstacle entrains and “holds-up” (traps) the contaminant material of the gas cloud. The hold-up mechanism within the wake of the obstacle is correctly predicted by the model. Generally, numerical predictions of concentration are seen to be in quite good agreement with the experimental data. It is noted that in interpreting the level of agreement between predictions and measurements, one needs to recognize that the measurements of the short time-averaged concentration obtained in trial numbers 26 and 29 correspond to one realization of the instantaneous dispersion, whereas the model predictions correspond to an ensemble-averaged concentration (and, as such, is associated with an average over an ensemble of realizations of the instantaneous dispersion). Because the model prediction is compared to one realization out of the ensemble of all possible realizations, the discrepancy between the model predictions and the experimental measurements cannot be smaller than the expected fluctuations of the individual realizations about the ensemble mean. In consequence, it is impossible for even a perfect model to give a greater precision in its predictions of the concentration than the expected concentration variability from realization to realization in the atmospheric dispersion. Keeping this caveat in mind, any remaining discrepancies between the predictions and the measurements could be probably attributed to the use of the standard k–ε turbulence model. It is well known that the standard k–ε model tends to over-predict levels of the turbulence kinetic energy for a flow impinging on a solid obstacle [the so-called the “stagnation point anomaly”, Durbin, 1996], resulting in prediction of a larger recirculation bubble in the wake region of a three-dimensional obstacle than is actually observed. For trial number 26, the more pronounced impingement of the approaching flow on the front face of the cylindrical gas tent resulting in a predicted larger recirculation bubble in the wake region of the tent will have the tendency to dilute the gas concentration faster after the tent collapses. In consequence, this will lead to a lower predicted concentration downstream of the cylindrical gas tent when compared to the associated measured concentration, as can be seen in Fig. 4. For trial number 29, a larger recirculation bubble behind the upwind cubical obstacle will entrain material from the dense gas cloud (that has spread upwind owing to the gravity slumping) into the wake earlier. This effect will result in an earlier predicted arrival time for the cloud (at the leeward face of the upwind obstacle) in comparison to the measurements, as can be seen in Fig. 5. The effect of the buoyancy production term Gk in the k–ε model [cf. Eqs. (17) and (18)] on the prediction of dense gas dispersion is investigated and demonstrated in Fig. 6. It was found that the inclusion of the Gk term in the k–ε model suppresses the turbulence kinetic energy, particularly in the wake region of the obstacle, which is not too surprising owing to the negative buoyancy effect of the dense gas cloud. In contrast, the effect of the Gk term on the prediction of the concentration levels in the cloud appears to be insignificant (at least for the example considered here). This observation is similar to the findings of Worthy et al. (2001), where it was reported that the Gk term made very little difference to the modeling of buoyant (hot) thermal plumes. 3.2. Pipeline rupture with axisymmetric hill(s) The experiment of Simpson et al. (2002) consists of measurements of the flow over a three-dimensional axisymmetric hill

K.-J. Hsieh et al. / International Journal of Greenhouse Gas Control 17 (2013) 127–139

135

Trial 29: t = 8.47 s

5

10

K 0.10 0.09 0.08 0.07 0.06 0.05 0.04 0.03 0.02 0.01

10

K 0.10 0.09 0.08 0.07 0.06 0.05 0.04 0.03 0.02 0.01

10

C 100.00 46.42 21.54 10.00 4.64 2.15 1.00 0.46 0.22 0.10

10

C 100.00 46.42 21.54 10.00 4.64 2.15 1.00 0.46 0.22 0.10

y/H

4 3 2 1 0

x/H

5

Trial 29: t = 8.47 s, Gk = 0

5

y/H

4 3 2 1 0

x/H

5

Trial 29: t = 8.47 s

5

y/H

4 3 2 1 0

x/H

5

Trial 29: t = 8.47 s, Gk = 0

5

y/H

4 3 2 1 0

x/H

5

Fig. 6. Buoyancy production term sensitivity analysis: contours of turbulence kinetic energy (top two panels) and mean concentration (% v/v) in the gas cloud (bottom two panels) in a vertical plane at z/H = 0 for trial number 29, with and without inclusion of the Gk term in the turbulence transport equations.

mounted on the floor of a wind tunnel. The wind tunnel has a working section that is 7.62 m long, 0.91 m wide, and 0.25 m high. The shape of the three-dimensional hill is defined as follows:





y(r) 1 r =− J0 ( )I0 H 6.04844 2H





− I0 ( )J0

r 2H



,

(20)

where H = 0.078 m is the hill height, = 3.1926, J0 (x) is the Bessel function of the first kind of order 0, and I0 (x) is the modified Bessel function of the first kind of order 0. The nominal inlet velocity for the spatial domain containing the axisymmetric hill is Uref = 27.5 m s−1 at the hill height H, and the Reynolds number based on H and Uref is 130,000. Measurements of the flow were conducted at various spanwise locations at 3.63H downstream of the hill crest. Detailed three-component laser Doppler anemometer measurements of velocity are available for this test case, providing a comprehensive documentation of a complex flow field with a challenging threedimensional separation. This experiment was used as “Test Case 11.2” in the 11th ERCOFTAC/IAHR Workshop on Refined Turbulence Modeling held at Chalmers University of Technology in Sweden on 7th and 8th April, 2005. The experimental data for this test case can be found at http://www.tfd.chalmers.se/∼gujo/WS11 2005/.

Simulations for flow over the axisymmetric hill were conducted on a computational mesh of 151 × 35 × 102 control volumes (nodes) in a spatial domain (shown in Fig. 7) with an extent of 19.8H × 3.205H × 11.67H in the streamwise, vertical and spanwise directions, respectively. Predicted results of the flow were compared with the associated experimental data at four spanwise locations (namely, at z/H = 0.00, 0.33, 0.65 and 0.81) at a downstream fetch of x/H = 3.63. At the inlet boundary of the spatial domain, the vertical profile of the mean streamwise velocity was approximated by a power-law



1/7

¯ function of the form u(y) = Uref y/H . At the outlet boundary, zero Neumann boundary conditions were imposed on all flow variables. At the spanwise boundaries, a symmetry boundary condition was used. Finally, at the top and bottom boundaries of the domain, standard wall functions were applied for mean velocities and turbulence quantities, and the zero Neumann boundary condition was applied for the mean concentration field. Fig. 8 shows the velocity vector field in the vicinity of the hill. In the experiment, there was a thin recirculation region observed on the leeward side of the hill. Simulations using RANS, large-eddy

136

K.-J. Hsieh et al. / International Journal of Greenhouse Gas Control 17 (2013) 127–139

y/H

6 4 2 0 -4

0

4

x/H

8

12

16 Y

-6 -4 -2

z/H 0

Z

X

2 4 6 -4 0 4

x/H

8 12

Fig. 7. A two-dimensional view of the computational mesh in a vertical x–y plane at z = 0 (top panel) and a three-dimensional view of the computational mesh along the ground surface (bottom panel). The mesh consists of 151 × 35 × 102 nodes in the streamwise (x), vertical (y), and spanwise (z) directions, respectively.

simulation (LES), or hybrid RANS-LES (e.g., Garcia-Villalba et al., 2009; Krajnovic, 2008; Persson et al., 2006; Tessicini et al., 2007) have shown that it is difficult to predict this recirculation region correctly. In the RANS calculation of Persson et al. (2006) in which a non-linear cubic eddy viscosity model was employed, these investigators predicted a recirculation region that was too large. In one of the LES studies conducted by Tessicini et al. (2007), no recirculation was captured when they used a dynamic subgrid stress (SGS) model on a computational grid of 192 × 64 × 128 control volumes for a spatial domain with an extent of 16H × 3.205H × 11.67H, with the near-wall nodes having values of the normalized wall-normal distance y+ ( yu /, where u is the friction velocity and  is the kinematic molecular viscosity of the fluid) in the range between 40 and 60. In conformance with these previous studies, the standard k–ε model with wall functions used herein also failed to predict the recirculation region in the lee of the hill.

3

y/H

2

1

0 -3

-2

-1

0

x/H

1

2

3

Fig. 8. Mean velocity vector field in a vertical x–y plane at z/H = 0 in the vicinity of the hill, where H is the height of the hill.

Figs. 9 and 10 show the predictions of mean streamwise velocity and turbulence kinetic energy, respectively, at various spanwise z-locations at a fixed downwind fetch of x/H = 3.63. In general, the predicted mean velocity profiles are in good conformance with the measurements. However, the predicted turbulence kinetic energy profiles in Fig. 10 show a poorer agreement with the measurements. For example, the measured turbulence kinetic energy at z/H = 0.0 decreases rapidly for y/H > 0.9, but the numerical prediction in this region is quite high. This is likely due to the standard k–ε model predicting excessive turbulence kinetic energy levels when the approaching flow impinges on the windward side of the hill. This high turbulence energy is advected downstream, affecting predictions of the levels turbulence kinetic energy at all locations downwind of the impingement zone. Although there is no concentration data available from Simpson et al. (2002) to compare with our predictions of dense gas dispersion on complex topography, the three-dimensional axisymmetric hill is used to study the effect of the terrain on the dispersion of a dense gas cloud. A scenario is proposed here to investigate the potential safety concern involving a pipeline rupture (through either an accidental or malicious release) in the vicinity of hilly terrain. It is assumed that CO2 is released continuously from a section of pipeline (on the ground surface) located at x/H = −3 upstream of the center (summit) of the hill. The area of the source release is 0.118H2 m2 and the gas jet velocity from the release is 0.1Uref m s−1 . Atmospheric dispersion of the resulting dense gas plume over one and two hills is considered in our simulations. Fig. 11 shows the computational mesh for the case of the two hills, where the distribution of grid nodes and domain size is identical to that for the one hill case. Fig. 12 shows the mean flow streamlines and isosurfaces of the concentration field resulting from the CO2 release. Here, two iso-surfaces of the concentration field with values of 0.5% and 4% (v/v) are selected to exhibit various features of the resulting

K.-J. Hsieh et al. / International Journal of Greenhouse Gas Control 17 (2013) 127–139

1.5 1.2

1.5

z/H=0.00

1.2 0.9

y/H

0.6

0.6

0.3 0 0

z/H=0.33

y/H

0.9

0.3 0.2

0.4

0.6

0.8

1

0 0

1.2

0.2

0.4

u/Uref 1.5 1.2

0.6

0.8

1

1.2

0.8

1

1.2

u/Uref 1.5

z/H=0.65

1.2

z/H=0.81

0.9

y/H

y/H

0.9

0.6

0.6

0.3 0 0

137

0.3

0.2

0.4

0.6

u/Uref

0.8

1

0 0

1.2

0.2

0.4

0.6

u/Uref

Fig. 9. Vertical profiles of the normalized mean streamwise velocity at four spanwise z-locations (namely, at z/H = 0.00, 0.33, 0.65 and 0.81) at a fixed downstream fetch of x/H = 3.63: () experimental data of Simpson et al. (2002); (—) simulation. Here, H is the height of the hill, and Uref is the streamwise mean velocity at hill height.

1.5 1.2 0.9 0.6

0.02

2 ref

k/U

1.5

0.6

0 0

0.06

0.9 0.6 0.3

0.02

k/U2ref

1.5

0.04

0.06

z/H=0.81

1.2

y/H

y/H

0.04

z/H=0.65

1.2

0 0

0.9

0.3

0.3 0 0

z/H=0.33

1.2

y/H

y/H

1.5

z/H=0.00

0.9 0.6 0.3

0.02

k/U2ref

0.04

0.06

0 0

0.02

k/U2ref

0.04

0.06

Fig. 10. Vertical profiles of the normalized turbulence kinetic energy at four spanwise z-locations (namely, at z/H = 0.00, 0.33, 0.65 and 0.81) at a fixed downwind fetch of x/H = 3.63: () experimental data of Simpson et al. (2002); (—) simulation. Here, H is the height of the hill, and Uref is the mean streamwise velocity at hill height.

138

K.-J. Hsieh et al. / International Journal of Greenhouse Gas Control 17 (2013) 127–139 Y

-6 -4 -2

z/H 0

Z

X

2 4 6 -4 0 4

x/H

8

Fig. 11. A three-dimensional view of the computational mesh along the ground surface for the case of two hills.

concentration distribution. As a point of reference, we note that the National Institute for Occupational Safety and Health (NIOSH) recommends chronic exposure limits for CO2 for a 40-h work week to not exceed 0.5% (5000 ppmv). Furthermore, NIOSH suggests that an acute exposure to CO2 concentration levels at 4% (40,000 ppmv) is immediately dangerous to life or health (NIOSH, 2005). The effect of terrain (viz., one hill compared to two hills) on the concentration field is displayed in Fig. 12 for concentration levels

of 0.5% and 4% (v/v). It can be seen that significant concentrations at a level of 4% (bright green contour) exist at or near the ground surface in the range of downwind locations lying between the locations of the pipeline rupture and the first hill (owing to the fact that the plume of dense gas is “pancake-shaped” and does not rise very far from the ground). For the case of one hill, it is seen that the concentration iso-surface of 0.5% (dark green contour) extends around and over the single hill topography to regions that are well downstream of the lee of the hill. In contrast, for the case of dispersion over two hills, the distribution of the concentration iso-surface at a value of 0.5% is more limited in its downwind extent. In particular, the presence of the second hill appears to “trap” the contaminant in the region (valley) between the two hills, resulting in an accumulation of CO2 here. Furthermore, there exists also a recirculation region on the leeward side of the second hill, and this region will facilitate further entrainment and holdup of the CO2 in the leeward region of the second hill.

4. Conclusions

Fig. 12. Mean flow streamlines along the ground surface and two iso-surfaces of concentration in the dense gas plume at concentration levels of 4% (bright green contour) and 0.5% (dark green contour) for dispersion over one hill (top panel) and two hills (bottom panel). The concentration is expressed in (% v/v). (For interpretation of the references to color in this figure legend, the reader is referred to the web version of the article.)

Most of the numerical studies on CCS have been focused on the subsurface modeling, with the consequence that only very few studies have investigated the leakage of CO2 from the subsurface into the atmosphere and subsequent dispersion of the dense gas cloud. In this paper, a dense gas dispersion model based on the RANS approach has been developed and applied to simulate two scenarios that serve as surrogates for potential CO2 leakage events from CCS related infrastructure, such as pipelines and CO2 storage tanks in complex terrain and/or in the vicinity of buildings and other obstacles and obstructions. A simple and rigorous derivation for the modeling of the buoyancy term within the framework of the Boussinesq approximation is provided. The effects of buildings and/or complex topography on the dispersion of dense gas clouds/plumes are studied. It is shown that the model is able to provide fair to good predictions for either the flow or species concentration for the two scenarios studied herein. The prediction of the distribution of the dense gas cloud concentration when expressed in terms of iso-contours or iso-surfaces of critical concentration levels provides useful information for the delineation of hazard zones (or, toxic corridors) resulting from the release of a dense gas cloud near ground level. It is anticipated that the model can be used to provide information required for undertaking quantified risk assessments for the potential risks to the public posed by the accidental or malicious release of carbon dioxide from existing or future CCS related infrastructure.

K.-J. Hsieh et al. / International Journal of Greenhouse Gas Control 17 (2013) 127–139

Acknowledgements This work was supported by CanmetENERGY – Ottawa, Natural Resources Canada. References Class, H., Ebigbo, A., Helmig, R., Dahle, H.K., Nordbotten, J.M., Celia, M.A., Audigane, P., Darcis, M., Ennis-King, J., Fan, Y., Flemisch, B., Gasda, S.E., Jin, M., Krug, S., Labregere, D., Beni, A.N., Pawar, R.J., Sbai, A., Thomas, S.G., Trenty, L., Wei, L., 2009. A benchmark study on problems related to CO2 storage in geological formations: summary and discussion of the results. Computational Geosciences 13, 409–434. Davies, M.E., Singh, S., 1985. The phase II trials: a data set on the effect of obstructions. Journal of Hazardous Materials 11, 301–323. de Lary, L., Loschetter, A., Bouc, O., Rohmer, J., Oldenburg, C.M., 2012. Assessing health impacts of CO2 leakage from a geological storage site into buildings: role of attenuation in the unsaturated zone and building foundation. International Journal of Greenhouse Gas Control 9, 322–333. Durbin, P.A., 1996. On the k-3 stagnation point anomaly. International Journal of Heat and Fluid Flow 17, 89–90. Ferziger, J.H., Peric, M., 2002. Computational Methods for Fluid Dynamics. SpringerVerlag. Garcia-Villalba, M., Li, N., Rodi, W., Leschziner, M.A., 2009. Large-eddy simulation of separated flow over a three-dimensional axisymmetric hill. Journal of Fluid Mechanics 627, 55–96. Hedlund, F.H., 2012. The extreme carbon dioxide outburst at the Menzengraben potash mine 7 July 1953. Safety Science 50, 537–553. Hosa, A., Esentia, M., Stewart, J., Haszeldine, S., 2011. Injection of CO2 into saline formations: benchmarking worldwide projects. Chemical Engineering Research and Design 89, 1855–1864. Koornneef, J., Spruijt, M., Molag, M., Ramirez, A., Turkenburg, W., Faaij, A., 2010. Quantitative risk assessment of CO2 transport by pipelines – a review of uncertainties and their impacts. Journal of Hazardous Materials 177, 12–27. Krajnovic, S., 2008. Large-eddy simulation of the flow over a three-dimensional hill. Flow, Turbulence and Combustion 81, 189–204. Lien, F.S., Leschziner, M.A., 1994. Upstream monotonic interpolation for scalar transport with application in complex turbulent flows. International Journal for Numerical Methods in Fluids 19, 527–548. Lien, F.S., Yee, E., Wang, B.C., Ji, H., 2010. Grand computational challenges for prediction of the turbulent wind flow and contaminant transport and dispersion in the complex urban environment. In: Lang, P.R., Lombargo, F.S. (Eds.), Atmospheric Turbulence, Meteorological Modeling and Aerodynamics. Nova Science Publishers. Mazzoldi, A., Hill, T., Colls, J.J., 2011. Assessing the risk for CO2 transportation within CCS projects, CFD modeling. International Journal of Greenhouse Gas Control 5, 816–825. Mazzoldi, A., Hill, T., Colls, J.J., 2008. CFD and Gaussian atmospheric dispersion models: a comparison for leakages within carbon dioxide transportation and storage facilities. Atmospheric Environment 42, 8046–8054. Mazzoldi, A., Picard, D., Sriram, P.G., Oldenburg, C.M., 2013. Simulation-based estimates of safety distances for pipeline transportation of carbon dioxide. Greenhouse Gases: Science and Technology 3, 66–83.

139

Michael, K., Golab, A., Shulakova, V., Ennis-King, J., Allinson, G., Sharma, S., Aiken, T., 2010. Geological storage of CO2 in saline aquifers – a review of the experience from existing storage operations. International Journal of Greenhouse Gas Control 4, 659–667. NIOSH, 2005. NIOSH Pocket Guide to Chemical Hazards, DHHS (NIOSH) Publication No. 2005-149. U.S. Government Printing Office. Ogretim, E., Gray, D.D., Bromhal, G.S., 2009. Effects of crosswind-topography interaction on the near-surface migration of a potential CO2 leak. Energy Procedia 1, 2341–2348. Oldenburg, C.M., Unger, A.J.A., 2004. Coupled vadose zone and atmospheric surfacelayer transport of carbon dioxide from geologic carbon sequestration sites. Vadose Zone Journal 3, 848–857. Patankar, S.V., Spalding, D.B., 1972. A calculation procedure for heat, mass and momentum transfer in three-dimensional parabolic flows. International Journal of Heat and Mass Transfer 15, 1787–1806. Perdikaris, G.A., Mayinger, F., 1993. Numerical simulation of heavy gas cloud dispersion within topographically complex terrain. Journal of Loss Prevention in the Process Industries 7, 391–396. Persson, T., Liefvendahl, M., Bensow, R.E., Fureby, C., 2006. Numerical investigation of the flow over an axisymmetric hill using LES, DES and RANS. Journal of Turbulence 7, 1–17. Pielke Sr., R.A., 2002. Mesoscale Meteorological Modeling, second edition. Academic Press. Pruess, K., Oldenburg, C., Moridis, G., 1999. TOUGH2 User’s Guide. Version 2.0, Report LBNL-43134. Lawrence Berkeley National Laboratory, Berkeley, CA. Pruess, K., Garcia, J., Kovscek, T., Oldenburg, C., Rutqvist, J., Steefel, C., Xu, T., 2004. Code intercomparison builds confidence in numerical simulation models for geologic disposal of CO2 . Energy 29, 1431–1444. Pruess, K., 2008. Leakage of CO2 from geological storage: role of secondary accumulation at shallow depth. International Journal of Greenhouse Gas Control 2, 37–46. Qin, T.X., Guo, Y.C., Lin, W.Y., 2007. Large eddy simulation of heavy gas dispersion around an obstacle. In: Proceedings of the Fifth International Conference on Fluid Mechanics, Shanghai, China. Rhie, C.M., Chow, W.L., 1983. Numerical study of the turbulent flow past an airfoil with trailing edge separation. AIAA Journal 21, 1525–1532. Rodi, W., 1978. Turbulence Models and Their Applications in Hydraulics – A State of the Art Review. University of Karlsruhe, SFB 80/T/127. Scargiali, F., Di Rienzo, E., Ciofalo, M., Grisafi, F., Brucato, A., 2005. Heavy gas dispersion modeling over a topographically complex mesoscale: a CFD based approach. Process Safety and Environmental Protection 83, 242–256. Schnaar, G., Digiulio, D.C., 2009. Computational modeling of the geologic sequestration of carbon dioxide. Vadose Zone Journal 8, 389–403. Simpson, R.L., Long, C.H., Byun, G., 2002. Study of vortical separation from an axisymmetric hill. International Journal of Heat and Fluid Flow 23, 582–591. Sklavounos, S., Rigas, F., 2004. Validation of turbulence models in heavy gas dispersion over obstacles. Journal of Hazardous Materials 108, 9–20. Tessicini, F., Li, N., Leschziner, M.A., 2007. Large-eddy simulation of threedimensional flow around a hill-shaped obstruction with a zonal near-wall approximation. International Journal of Heat and Fluid Flow 28, 894–908. Worthy, J., Sanderson, V., Rubini, P., 2001. Comparison of modified k-ε turbulence models for buoyant plumes. Numerical Heat Transfer, Part B 39, 151–165. Yee, E., Lien, F.S., Ji, H., 2007. Technical description of urban microscale modeling system: component 1 of CRTI Project 02-0093RD. DRDC, Suffield, TR 2007-067.