Simulationbased estimates of safety distances for pipeline ...

See discussions, stats, and author profiles for this publication at: https://www.researchgate.net/publication/280041167

Simulation-based estimates of safety distances for pipeline transportation of carbon dioxide Data · July 2015

CITATIONS

READS

3

29

1 author: Alberto Mazzoldi ----20 PUBLICATIONS 287 CITATIONS SEE PROFILE

Some of the authors of this publication are also working on these related projects:

Potential fields analysis for geothermal exploration View project

CCS Certification Framework View project

All content following this page was uploaded by Alberto Mazzoldi on 14 July 2015. The user has requested enhancement of the downloaded file. All in-text references underlined in blue are added to the original document and are linked to publications on ResearchGate, letting you access and read them immediately.

Modeling and Analysis

Simulation-based estimates of safety distances for pipeline transportation of carbon dioxide Alberto Mazzoldi, Earth Sciences Division, Lawrence Berkeley National Laboratory, Berkeley, CA David Picard, Papagudi G. Sriram, Clearstone Engineering Ltd., Calgary, Alberta Curtis M. Oldenburg, Earth Sciences Division, Lawrence Berkeley National Laboratory, Berkeley, CA Abstract: Pipeline transportation of fluids is a proven technology for moving large quantities of liquids and gases (e.g. hydrocarbons, hazardous liquids, hydrogen). The anticipated introduction of largescale geologic carbon sequestration (GCS) as a means of reducing greenhouse gas (GHG) emissions will require the ability to transport massive amounts of carbon dioxide (CO2) safely and economically. To accommodate GCS demands, the existing US and European CO2 pipeline infrastructure may eventually have to be expanded to be comparable in size to natural gas and oil pipeline systems. Furthermore, the new pipeline infrastructure will inevitably intersect with population centers as it connects sources with sink areas. There are important unanswered questions about pipeline network requirements, regulations, utility cost recovery, economics, regulatory classification of CO2 itself, and pipeline safety. The focus of this research is on this last aspect, i.e. safety of the general public, workers, and property related to the transportation of CO2. We carried out simulations that coupled two computational fluid dynamics (CFD) codes to determine: (i) leakage rates from fully ruptured aboveground CO2 pipelines for a typical pipeline fluid composition, and (ii) the resulting atmospheric dispersion of the gas near the broken pipe. Using threshold values for atmospheric CO2 concentration, our work shows that concentrations dangerous to human health can extend on the order of hundreds of meters from the ruptured pipeline. This work contributes to the knowledge base needed to establish safety distances for routing CO2 pipelines through inhabited and other sensitive areas. © 2012 Society of Chemical Industry and John Wiley & Sons, Ltd Keywords: geologic carbon sequestration; CO2 transportation; pipeline leakage risk; pipeline decompression; atmospheric dispersion; safety Cμ

Nomenclature a Cε1 Cε2 Cε3

s–1]

velocity of sound in a fluid [m k-ε turbulence model dimensionless constant equal to 1.44 k-ε turbulence model dimensionless constant equal to 1.92 k-ε turbulence model dimensionless constant equal to –0.1

Cf CP Dm Dp

dimensionless turbulence viscosity constant for the k-ε model, equal to 0.09 Fanning pipeline friction factor, dimensionless specific heat at constant pressure [J kg–1 s–1] effective diff usion coefficient for species m [m–2 s–1] pipe diameter [m, inches]

Correspondence to: Alberto Mazzoldi, Earth Sciences Division, Lawrence Berkeley National Laboratory, 1 Cyclotron, Road Berkeley, CA 94720 USA. E-mail: [email protected] Received June 11, 2012; revised October 3, 2012; accepted October 4, 2012 Published online at Wiley Online Library (wileyonlinelibrary.com). DOI: 10.1002/ghg.1318

© 2012 Society of Chemical Industry and John Wiley & Sons, Ltd | Greenhouse Gas Sci Technol. 2:1–18 (2012); DOI: 10.1002/ghg

1

A Mazzoldi et al.

g hm L k kc MWm P Pk Pb q R ℜ Sm Sp SU ST Skv t T T∞ U u u* ε

ζ φ к κ ym σ σk σε μ μt θ* ∇

gravitational acceleration, 9.81 m s–1 specific enthalpy of species m [J kg–1] Monin-Obukhov turbulent length scale [m] turbulent kinetic energy per unit mass [m2 s–2] effective thermal conductivity [W m–1 K–1] molecular weight of species m [kg] pressure [atm, MPa] mechanical production rate of k [J s–1] buoyancy production rate of k [J s–1] heat flux vector [W m–2] mixture of gas constant [J/mol–1 K–1] universal gas constant, 8.314 J/mol–1 K–1 source term for species m [kg m–3 s–1] source term for continuity [kg m–3 s–1] source term for momentum [kg m–3 s–1] source term for temperature [kg m–3 s–1] source term due to vegetation canopy and interaction with particles [kg m–3 s–1] time [s] temperature [K] temperature around the pipe [K] velocity vector [m s–1] fluid velocity [m s–1] surface friction velocity [m s–1] dissipation of turbulent kinetic energy per unit mass [m2 s–3] Monin-Obukhov similarity variable = z/L, dimensionless Joule-Thomson coefficient, φCO2 = 13 K MPa–1 dilatational viscosity [m2 s–1] Von Karman constant = 0.41, dimensionless mass fraction of species m Newtonian viscous stress tensor [N m–2] Prandtl number for turbulent diff usion, equal to 1.0 dimensionless turbulence model constant for the ε equation, equal to 1.314 effective viscosity [N s m–2] eddy viscosity [N s m–2] temperature scale nabla operator

Introduction t is widely believed that geological carbon sequestration (GCS) will be necessary to substantially reduce global CO2 emissions.1 GCS is also fundamental to achieving CO2 emissions reductions at the

I 2

Modeling and Analysis: Pipeline transportation of CO2

lowest cost: without it, the overall cost for reaching the required emission reductions could rise by up to 70%.2 Coal-, oil- or gas-fired power plants are the most likely initial candidates for using this technology because they constitute large single-point sources of CO2 and they contribute a large share of global emissions from fossil fuels.3 Today’s carbon capture technologies can potentially capture 80–95% of CO2 emitted from an electric power plant or other industrial source. Captured CO2 would then be transported to places suitable for underground injection via high-pressure pipelines. Approximately 1.5 billion tons of CO2 are produced annually in the USA from coal-fired power plants.4 If all of this CO2 were to be transported for sequestration, the quantity would be equivalent to three times the mass and, under typical operating conditions, one-third the volume of natural gas transported annually by the US natural gas pipeline system.5 Transporting the quantity of CO2 required by a widespread implementation of GCS would likely require a dedicated pipeline network.6 Pipeline transportation of CO2 at high pressure (≥8 MPa) is not without risk,7 the main hazard being represented by the sudden rupture of a pipeline that may lead to potential asphyxiation of nearby people and animals, CO2 being lethal if inhaled for approximately 10 min at concentrations larger than 105 ppm in air (10% by volume).8,9 Although in this study we focus on CO2 as a potentially dangerous gas, another relevant hazard associated with high-pressure pipeline transportation is the shock wave that forms by the expansion of gas as it exits from high-pressure containment and which is capable of killing people in the immediate proximity of the rupture, depending on the diameter and pressure of the pipeline.10 Based on models of costs for electricity transmission lines, CO2 pipelines, and fuel transportation systems, it is always preferable to locate a power plant with CO2 capture capability near the electric load (i.e. electricity end-users), reducing the losses and costs of bulk electricity transmission.11 This means CO2 pipelines will likely be routed near populated areas. CO2 pipeline leakage hazard is defined as the probability of occurrence of dangerous CO2 concentrations per year at any given distance from the hazard source. It is usually expressed in the form of iso-hazard contours around the hazard source and, in the case of pipelines, iso-contours are usually parallel to the pipeline route.12 Following the Precautionary



Principle, we consider individual full-bore ruptures of above-ground pipelines, and we evaluate the maximum distance away from ruptured pipelines that can be reached by hazardous CO2 concentrations. This is done within the limitation of not considering the under-expansion characteristic of jet flow after exiting high-pressure containment,13 which results in an over-prediction of CO2/air mixing near the source. This neglect of under-expansion leads to under-prediction of CO2 concentrations farther from the source, notwithstanding the stable atmospheric conditions, constant wind velocities, and lack of modeling of dry-ice formation (solid CO2) all of which tend to over-predict CO2 concentrations downwind from the pipe rupture as we will discuss later. During the last two decades, several studies considered atmospheric dispersion of CO2 leaking from GCS transportation facilities.14–17 These studies map hazard using maximum downwind extents of hazardous concentrations of CO2, evaluated through the use of Gaussian/dense-gas models to simulate dispersion of gas plumes in the atmosphere, assuming constant release rates during the decompression of the pipeline. These previous studies did not take into account the initial momentum of the gas as it exits the damaged pipeline, or the time-varying flow-rates of leakage from the rupture. Results from previous research highlight how hazard evaluation through the usage of Gaussian models can be overly conservative (lead to unrealistically large hazard areas).18 Computational fluid dynamics (CFD) atmospheric dispersion models numerically integrate the conservation equations for mass, momentum, and energy of fluid flows in order to evaluate pressure, temperature, density, and energy of a gas dispersing in time and space. CFD codes are the most rigorous treatments of the physics of fluid flow and have been used for simulating the high-speed release of flammable gases like natural gas (CH4) and hydrogen due to leaks from high-pressure transportation facilities.19–22 One of the key controls on the ultimate atmospheric dispersion of leaking CO2 is the source term. Fundamental parameters characterizing the high-pressure CO2 release source term and thereby influencing the dispersion of the gas in the atmosphere are: (i) the composition of the gas; (ii) the mass-flow rate and velocity of the gas exiting the pipe (at the source), and (iii) the thermodynamic variables P and T. In this study, we use the CFD non-isentropic, pipeline-blowdown PIPE model,23 to characterize

A Mazzoldi et al.

transient CO2 releases after sudden full-bore rupture from pipelines of different dimensions (diameter and length) at Pi = 100 bar and transporting a fluid mixture of 97% CO2, 2% CH4 and 1% N2. The chosen composition is a simple three-component mixture representative of CO2 injection mixtures encountered in tertiary enhanced oil recovery (EOR) operations involving CO2 flooding. We simulate the process of leakage from one end of a ruptured pipe and associated depressurization of the pipe. The scenario of a full-bore rupture of a pipeline would, by definition, result in two broken-end pipe segments releasing fluid in different directions. For simplicity in this study we consider only one of these broken ends, a scenario that would be representative of a pipeline break near an emergency shutdown (ESD) valve that would interrupt leakage from one of the broken ends immediately after the rupture but not from the other. Our approach can be easily extended to simulate releases from two (or more) broken-end arrangements. We use pipe decompression data obtained from the PIPE model as input to the CFD atmospheric dispersion code Fluidyn-PANACHE which has been applied previously for atmospheric dispersion related to pipeline risk assessment.18,24 We investigate various pipeline dimensions and how they affect safety issues for CO2-GCS transportation.25 In this study, we give estimates of maximum distances covered by dangerous concentrations of the dense gas (defined as gas that is heavier than air at any P-T conditions). For conservatism that maximizes safety distances, we consider horizontal releases with a direction perpendicular to the pipeline. Results of this research can be used in determining the routes of new pipelines and safety issues involved in its planning (e.g. diameter of the pipe, ESD locations and frequency, distance from highly sensitive locations, and ground cover).

Hazard and risk assessment Risk is commonly defined as the product of the likelihood of a hazardous event and its consequences:26 Risk = Hazard × Consequences

(1)

The main hazard involved in the transportation of CO2 is accidental leakage from transportation facilities. Consequences are the anticipated effects of high


3

A Mazzoldi et al.


Table 1. Yearly failure rate summary per module.17 Module

1

CO2 recovery at source

Failure rate (per module per year)

Leak every x years, x=

1.5 * 10–1

7

–3

2

Converging pipelines

4.6 * 10

217

3

Booster station

4.0 * 10–2

25

10–4

4

Pipelines

3.4 *

5

Injection well

1.8 * 10–1

2,941 6

CO2 concentrations in air on human health and the environment. Carbon dioxide is toxic to humans and other animals at high concentrations. The Immediately Dangerous to Life and Health (IDLH) limit set by the National Institute for Occupational Safety and Health (NIOSH) is 40 000 ppm (4% by volume). At concentrations larger than 10% by volume (100 000 ppm), CO2 may cause unconsciousness after a few minutes and would be potentially fatal after 10–15 min of continuous inhalation.17 At concentrations above 25% (250 000 ppm) it poses a significant risk of asphyxiation, impeding motion after a few breaths and killing after less than 1 min of exposure.27 Transportation of CO2 by pipeline has been carried out since the early 1970s for EOR projects.28 The existing database for CO2 is relatively small and pipeline failure frequencies are borrowed from the hydrocarbon transportation systems experience.8,17 Table 1 presents assumed overall failure frequencies for the different modules of the CO2 transportation systems, considering a range of potential leakage dimensions. Consequences of a hazardous event are difficult to estimate and generalize. The effects on humans from a given release of CO2 in an area depend primarily on the total amount of CO2 leaked, its concentration in the atmosphere, and the population density in the area.* A thorough hazard analysis would cover the whole range of hazardous events (leak dimensions) to obtain a model for the total hazard of pipeline failure. Here, we limit our study to the analysis of a conservative scenario, consistent with application of the Precautionary Principle.10,29 By this means, in order to achieve results that can be of help in the determi* CO

2 being denser than air, its concentration in the atmosphere would also depend on the topography of the area. This parameter is site-specific and, in order to achieve general results, we have limited our study to flat and horizontal terrain.

4

nation of safety distances from CO2 pipelines, we focus our attention on full-bore ruptures, which we consider worst-case scenarios. We then determine the downstream safety length (DSL) as the maximum distance away from the source (pipe rupture) contained by the two dangerous iso-concentration lines, i.e. 100 000 and 250 000 ppm.

Special considerations for high pressure CO2 releases There are four different phases for CO2: solid, liquid, gas, and supercritical (Fig. 1). Above the critical pressure and temperature (P = 7.38 MPa and T = 31 °C), CO2 is a supercritical fluid with a liquidlike density that also behaves like a gas in occupying the volume of its container.30 Typical CO2 pipelines operate at temperatures ranging between 12 °C and 44 °C and pressures ranging between 85 bars and 150 bars.31 CO2 is normally transported in the liquid phase,28 although it can transition without significant change in volume or other properties to its supercritical form if the ambient (pipeline) temperature slightly exceeds 31 °C. When leaking from high-pressure containment, the large pressure differential experienced by CO2 will cause it to eventually transition from liquid to gas, with associated density decrease and consequent high velocity when exiting the damaged pipeline. The compressibility factor (Z = PV/nRT) of CO2 is highly non-linear and, at P = 10 MPa, Z < 0.3.32 During decompression, the Joule-Thompson effect will drastically reduce the temperature according to: ΔT = φ · ΔP

(2)

where ϕ is the Joule-Thomson coefficient. For CO2, the value of ϕ is ϕCO2 = 13 K MPa–1.32 For the cases considered in this paper, with an initial P = 10 MPa, the maximum temperature drop caused by the pressure change (~10 MPa) would be around 130 K: CO2 molecules may experience a phase change to solid state (dry ice; Figs 1 and 2(d)).33 At near-sublimation temperature (–78 oC) gaseous CO2 has a density of about 2.8 kg m–3, higher than at ambient P-T conditions (ρCO2 = 1.8 kg m–3) and much higher than air (ρair = 1.2 kg m–3). Figure 2 shows how the high-momentum-flow of gas into a static atmosphere encounters resistance and forces air entrainment into the plume. Frictional heating plus the ambient air temperature provide heat



A Mazzoldi et al.

Figure 1. CO2 phase diagram.30 1 atm = 0.1013 MPa.

to the CO2, which can convert from solid back to gas. Upward and horizontal releases do not involve large losses of mass from the flow in the form of dry-ice deposition on the ground (as happens for downward releases, Fig. 2(d)).33 After losing its momentum due to air-resistance, CO2 continues to mix with air, facilitated by wind and atmospheric turbulence.34

CFD models and simulations For modeling the temporal evolution of hazardous concentrations of CO2 over an area after the start of the release, the CFD atmospheric dispersion model PANACHE is fed with source-term data (leakage rate, flow velocity, and temperature of the jet-flow, as they vary in time after beginning of the release) derived from the PIPE simulations, considering pipelines of different dimensions (diameter and length) transporting a non-pure CO2 stream representing an injection mixture used in CO2-EOR activity. Below, we describe the CFD models, and give a summary of the main assumptions and results of maximum distances away from pipes reached by dangerous concentrations in the atmosphere of the gas in question (10% and 25% by volume).

PIPE decompression model Transient leakage rates from high-pressure CO2 pipelines were predicted using the one-dimensional,

non-isentropic, real-fluid blow-down simulation code developed by Picard and Bishnoi.23 This code is based on the one-dimensional transient Navier-Stokes mass, momentum, and energy balance equations uniquely expressed in terms of fluid velocity, density, and pressure to facilitate convenient implementation of a ‘real fluid’ equation of state: D ρ ∂u Mass Conservation Equation: (3) + =0 Dt ∂xx Momentum Conservation Equation:

ρ

μ|μ| Du ∂P −2C f ρ + = Dt ∂xx Dp

(4)

Energy Conservation Equation:

Dp Dρ e − a2 = Dt Dt

C f ρu3 + 4U (T° − T )

D p ρT

(5)

In these equations, t and x denote, respectively, time since the start of the leak and axial position in the pipeline, Dp denotes the pipe diameter and ρ, P, T and a are local fluid density, pressure, temperature and sound velocity, respectively. The term D/Dt is the substantial differential operator ∂/∂t + u ∂/∂x. Cf is the Fanning friction factor.39 T∞ is the temperature around the pipe. The term e denotes the term ⎛⎜ ∂p ⎞⎟ . ⎝ ∂s ⎠ ρ


5

A Mazzoldi et al.


Figure 2. a) A schematic representation of how the jet-mixing effect operates on a high momentum flow of CO2 from a leaking pipeline;35 b) an upward jet release from a HP surface facility, in which air provides enough heat for CO2 to reconvert into gas (Barendrecht, NL);36 c) a horizontal jet-flow of CO2 is not expected to lose much material in the form of dry ice (Cranfield, MS);37 d) a downward release would give rise to the accumulation of dry ice on the ground (Sutton-Bonnington Campus, University of Nottingham, UK).38

The solution of these equations assumes the following: (i) the fluid flow is one-dimensional; (ii) frictional forces are acting only on the fluid on the inside surface of the pipe wall; (iii) there is heat exchange through the pipe wall with the surrounding

6

environment; and (iv) the fluid is always in thermodynamic equilibrium. For improved numerical accuracy, the equations are transformed, using the Method of Characteristics, into ordinary differential equations that were solved



using the implicit finite-difference scheme described in Chaudhry.40 The equation of state used to estimate the thermodynamic properties and phase behavior of CO2 is the multi-component equation of state (EOS) of Peng and Robinson (1976),41 which includes the capability of modeling impurities in the CO2 (such as N2, H2S, and CH4).

The Peng-Robinson equation of state To justify the use of the Peng-Robinson EOS for pipeline blowdown involving potentially very low temperatures, we offer the following. Figure 3 presents a comparison of the decompression curves for the three-component CO2 mixture as modeled using the two-phase (i.e. vapor and liquid) Peng-Robinson EOS to the corresponding curves determined for pure CO2 using the three-phase (i.e. vapor, liquid and solid) Span and Wagner EOS.42 The Span and Wagner EOS predicts the solid-phase thermodynamic properties using a sublimation curve and a solid CO2 specific heat of 1.1 kJ/kg/K,43 and density of 1562 kg/m3.44 There is no significant difference in temperature predictions shown in Fig. 3(a) between the two thermodynamic models despite the presence of minor amounts of impurities in the three-component CO2 mixture used with the Peng-Robinson EOS. In Fig. 3(b), the difference in critical mass flux (mass flux at Mach 1 flow) is not significant in the two-phase region. However, there is significant difference between the two models in the supercritical region. The prediction of density using the Peng-Robinson EOS is lower than that for pure CO2 as shown in Fig. 3(c). Figure 3(d) shows the predicted sound speed during decompression using the two equations of state. Figure 3(e) presents the predicted quality (vapor-phase mass fraction) of the CO2 stream during the decompression process. In the supercritical region, there is only one phase and that phase is considered to be vapor-like when the derivative of thermodynamic sound speed with respect to temperature is greater than or equal to zero, and, it is liquid-like when it is less than zero. The Peng-Robinson equation of state for the three-component system predicts a uniform liquid-like characteristic in the supercritical region and increasing vapor formation with reduction in pressure in the two-phase region. The pure CO2 EOS predicts vapor-like behavior initially and as the pressure approaches the critical pressure of pure CO2

A Mazzoldi et al.

the supercritical phase changes to liquid-like behavior. In the two-phase region the pure CO2 equation of state initially predicts a vapor-liquid mixture, but the liquid phase eventually turns to a solid phase as the pressure decreased below the triple-point pressure, while the Peng-Robinson EOS predicts a vapor-liquid mixture throughout and a greater vapor phase mass fraction when the pressure is above the triple point pressure of pure CO2. At the triple point pressure (i.e. 517 kPa), the pure CO2 decompression curve for quality shows a horizontal discontinuity in the predicted vapor phase mass fraction and the predicted vapor-phase mass fraction is greater than the predictions with the Peng-Robinson EOS (Fig. 3(e). For the two-phase equilibrium conditions below the triple point pressure of CO2, the Peng-Robinson EOS predicts a lower amount of vapor when compared to the pure CO2 EOS and does not show any discontinuities in this same region. Despite the differences in density and sound speed predictions in the two-phase region and the fact that the Peng-Robinson EOS cannot predict the presence of a solid phase, the critical mass flux prediction is not significantly different from that predicted using the pure CO2 EOS. This indicates that the mass flow rate predicted using the Peng-Robinson EOS is reliable for choked flow in two-phase conditions (i.e. the effects of overestimating the sound speed and underestimating density tend to cancel out to result in reasonable mass flow rate predictions), which characterizes most of the decompression process.

PANACHE atmospheric dispersion model Fluidyn-PANACHE (version 4.0724),45 is a computer code for numerical simulation of 3D atmospheric flows and pollutants dispersion in short and mediumrange scales. It solves the Navier-Stokes Equations (NSE) for the motion of the fluid and advection-diff usion transport equations for the species in the Reynolds Average Navier-Stokes (RANS) frame. The mean flow, temperature, and the species concentration are calculated from the RANS equations. From the non-linear terms in NSE, second order moments (e.g. Reynolds stresses) are formed, involving turbulent fluctuations (sub-grid scales and short term) and are responsible for momentum, heat, and mass turbulent fluxes. In Fluidyn-PANACHE, a linear eddy viscosity model (and heat and mass diff usivities model) is used as the closure scheme for such fluxes.


7

A Mazzoldi et al.


Figure 3. The predicted change in selected thermodynamic properties of pure carbon dioxide and a three-component carbon dioxide mixture during the decompression of a 15 km segment of 6 NPS pipeline from an initial pressure of 100 bars and temperature of 37 °C.

8



A Mazzoldi et al.

The related turbulent parameters are estimated from the average statistics of the turbulent fields in k, the turbulent kinetic energy and ε, the turbulent energy dissipation rate. Both mean flow equations and turbulent closure equations are described.

Mean flow and transport equations PANACHE solves the Navier-Stokes equations along with the equations describing conservation of species concentration, mass, and energy for a mixture of ideal gases. PANACHE solves the Reynolds averaged forms of these equations for turbulent flow. The Reynolds stresses are modelled using the linear eddy viscosity model.46 The governing equations are written below.47 Conservation of species concentration: ∂( ρ ym ) + ∇ ⋅ ( ρ ∪ ym ) = ∇ ⋅ Dm ∇( ym ) + Sm ∂t

(6)

where m = 1, n Summing equation (6) over all species, we get the continuity equations, or, ∂ρ Conservation of mass: (7) + ∇ ⋅ ( ρ ∪ ) = Sρ ∂t Conservation of momentum (N-S eq.):

∂

(ρ ) + ∇ ⋅ ( ρ ∪ ∪) + ∇ ⋅ τ − ∇P + S∪ ∂t

(8)

Conservation of energy: ⎡ ∂( ρT ) ⎤ Cp ⎢ + ∇ ⋅ ( ρ ∪ T )⎥ = ⎣ ∂t ⎦ ⎡ ∂P ⎤ − ∇ ⋅ q = ⎢ + ∪ ⋅ ∇P ⎥ + ⎣ ∂t ⎦

equation (same sum as for Sm), SU source term for momentum equation (sum of momentum flux related to the sources + drag force due to vegetation canopy + drag force on particles + body forces (e.g. gravitational acceleration)), ST source term for temperature equation (sum of total energy flux from sources + heat transfer to particles/droplets + energy release due to chemical reactions + heat source on boundaries (e.g. ground, water surfaces, domain borders)). Added to the mean flow and mean transport equations, the relevant thermodynamic equation models are used, particularly, ideal gas law is used for the thermodynamic model of mixture of gases and expression for specific heat capacity, Cp(T): P = ρRT

where, R = ℜ∑m ym/MWm mixture gas constant, ℜ = universal gas constant (8314.34 J/kmol.K), MWm = molecular weight of species m.

Turbulence model, domain size, and boundary conditions The standard k-ε model for the sub-grid turbulent closure for the RANV equations has been used throughout the simulations.48,49 Turbulent kinetic energy (k) and its dissipation rate (ε) are used in the turbulent viscosity and diff usivity (for mass and heat) operating in the mean flow equation:

μt =

∇ ∪ + ST

(9)

where, ρ denotes density (kg m–3), ∪ is velocity vector, ym mass fraction of species m, T is temperature (oC), P represents pressure, Dm effective diff usion coefficient . for species m, τ, equals to μ γ – (2/3μ – κ) (∇·U)δ is . the viscous stress tensor where γ = ∇U+(∇U)T is the rate-of-strain (or rate-of-deformation) tensor, μ is the effective viscosity, к denotes dilatational viscosity (= 0 for a Stokesian fluid), δ is unit tensor, Cp is specific heat at constant pressure, q represents the heat flux vector = –kc∇T, where kc is the effective thermal conductivity, Sm is the source term for species m equation (sum of pollution source + dry deposition on vegetation), Sρ is the source term for continuity

(10)

ρCμ K 2

(11)

ε

where Cμ is a constant equal to 0.09. The equations for turbulence generation and dissipation are given by: ⎛ μ ⎞ ∂( ρk ) + ∇ ⋅ ( ρ ∪ k ) = ∇ ⋅ μ + t ⎟ ∇k ∂t σk ⎠ ⎝ (12)

2 + Pk + Pb − ρε − ρk∇ ⋅∪ Skv 3 ⎛ μ ⎞ ∂( ρ ) ε + ∇ ⋅ ( ρ ∪ ) = ∇ ⋅ μ + t ⎟ ∇ε + [ ∂t σε ⎠ k ⎝ − Cz 2

]

2 C − Cε 3 ) ρε∇ ∪ Sεv 3 ε1

ε1

( Pk + Pb ) (13)

. where: Pk = μt γ : ∇∪ is the mechanical production rate of k, Pb = μtβ(g∇T/σk) is the buoyancy production rate of k, σk = 1.0 is the Prandtl number for turbulent diff usion of k, β is a constant equal to 20.0, Skv is the


9

A Mazzoldi et al.


source term due to vegetation canopy and interaction with particles. μt is computed using equation (11). Cs1, Cz2 Cs3 and σs are model constant whose value is given in the Nomenclature. The k- ε model computes the length scales from the local turbulence characteristics. Thus it can model flows subjected to both mechanical shear (obstacles, terrain undulations, canopy) as well as buoyancy (stability and buoyant/heavy gas plumes). The k-ε model is an isotropic model of turbulence that results in turbulent diff usivities that are same in both horizontal and vertical directions at a location, which may overestimate the effects of turbulence on the dispersion of a dense gas, but only in the far-field after the flow has lost its momentum. The domain sizes used throughout this study depend on the release rate considered for each numerical experiment (see below). Scenarios (domain size and mesh created) were linked to pipeline diameter: the three grids considered for the three scenarios (pipeline diameter (inches) of 6 in, 16 in, and 32) in all had between 30 000 and 80 000 cells, covering areas of 0.12 to 2 km2, with domain heights of 200 m. PANACHE can use a polygonal mesh generator that results in a finer grid where simulations need to be more precise (e.g. near a high-speed source). In simulations for this study, a nested domain has been used with a much finer grid around the source and downstream of it for all the scenarios considered, in order to capture the mixing of high-velocity CO2 and static air (in Fig. 6, the 123 m scale-segment represents the length of the nested domain downstream for the particular trial), which clearly resulted in a coarser mesh further away from the source. Properties on the surfaces of the domains are specified by boundary conditions. Ambient mean wind speed and air temperature profiles are specified within the model domain and are represented by logarithmic functions with corrections for the Atmospheric Boundary Layer (see PANACHE manual for further specifications),45 such that: v(z) = u*/κ[ln(z/z0) Ψ0(ζ)]

(14)

θ(z) = θ*/κ[ln(z/z0) Ψ2(ζ)]

(15)

where θ* = temperature scale; κ is the Von Karman constant equal to 0.41; Ψ1(ζ) and Ψ2(ζ) represent similarity profile. The surface friction velocity, u*, the temperature scale θ*, and the Monin-Obukhov length,

10

L are related by: L = u*2T / (g κ θ*) and θ* = Qh / (ρ Cp u*). The micrometeorological parameters u*, θ*, and L are evaluated for different atmospheric stability classes (see PANACHE User manual).45 In the simulations for this study, a wind velocity of 2 m s–1 at 10 m height and a T = 5 oC at ground level were used. The ground roughness is another boundary condition that can vary greatly with the nature of the area considered (presence of fields, forests, water bodies, etc.) and the presence of relatively small protruding objects. In this study, PANACHE was given a surface roughness length z0 = 0.001, typical of a lawn. After the creation and refining of the mesh, a typical simulation begins with the resolution of the wind field. Releases start a certain time after the beginning of the simulations (100s of seconds), when all the residuals of converging parameter (pressure, velocities, turbulence creation, etc.) are below a threshold, which for our study was set at 10 –2. During the first seconds after the start of a release, residuals increase, recording unbalance of the parameters, after which they typically settle down below 10–2.

CFD modeling assumptions In line with the Precautionary Principle, we modeled CO2 leakage from full-bore ruptures only. As for pressurized transport of other gases (CH4, H2), decompression cooling can result in increased embrittlement and sudden pipe failure, which can turn small leaks into large/full-bore ruptures. In the decompression simulations using PIPE, pipelines were considered as transporting a non-pure gas, the flow being made up of 97% CO2, 2% CH4 and 1% N2. According to the Peng-Robinson EOS, the above mixture at P = 100 atm and T = 15 oC, as in the simulations performed, has a density ρ = 827.5 kg m–3 and is liquid. Outputs of CFD simulations with the PIPE model were used as input to PANACHE simulations (see next paragraph). High-velocity source terms were modeled as horizontal velocities, at a height of 2 m above the ground surface. Only a single pipe release was simulated, implying that the leakage is occurring close to an ESD valve which shuts down the flow from the matching full-bore opening, immediately after the rupture. When leaking from a high pressure pipe, an underexpanded jet would form,13 a feature which we did not explicitly model. Neglect of under-expansion in atmospheric dispersion simulations can lead to over-prediction of CO2 dispersion in the near field and



consequent under-prediction of distances covered by specific CO2 concentration (e.g. 100 000 and 250 000 ppm, as in this work). In this study, we assume gas expansion occurs immediately after the gas exits the pipeline and we specify jet-flow rates for CO2 at the source (pipeline rupture) at ambient pressure. Furthermore, we consider that full-bore releases from buried lines would uncover the pipe soon after the start of the release,16 obviating the need to model leakage from underground pipes. In order to give a conservative evaluation of safety distance from a given pipeline, we modeled the leakage as taking place perpendicularly to the direction of the pipeline (see jet-flow velocity vectors in Fig. 6(a)). In terms of atmospheric conditions, it is known that atmospheric dispersion will increase with wind and a less stable atmosphere. Within every single run, the wind has been chosen with a direction parallel to the jet flow and at a low speed of 2 m s–1 (approximately 4.4 mph) at a 10 m height, considering a logarithmic wind profile (Fig. 4). The Pasquill-Gifford atmospheric stability was specified as class F, very stable,50,51 for all runs. We also modeled two experiments with wind having an opposite direction relative to gas release (Fig. 6). With the aim of leaving results open and applicable to different CO2 transportation projects, no particular geographic location was considered.

Source term characteristics The different pipeline geometries considered in the decompression simulations accounted for every combination of pipe-segment lengths and diameters, with values of 0.5, 1, 5, 10, 20, and 40 km and 0.15, 0.4, and 0.8 m (6, 16 and 32 inches), respectively. These values were considered in order to account for

Figure 4. Logarithmic wind profile

A Mazzoldi et al.

transportation pipelines over different distances and for diverse CO2 mass flow rates. Figure 5 shows results from PIPE simulations of leaking CO2 in time after rupture on four characteristic pipeline segments with different dimensions. The figure also shows our methodology for feeding source-term data from PIPE to PANACHE. In the figure, the blue lines represent the release rates as simulated by PIPE, while the red lines, superimposed in Figs 5(a) and 5(b), exemplify release rates used as input to PANACHE. Discrete rate values were used in PANACHE in order to decrease considerably the simulation times. The time length of each step within each release (specifically for the four first steps after beginning of the release) was one and the same for all of the simulations (i.e. t0 = 0, t1 = 2, t2 = 10, t3 = 30, and t4 = 80 s). Release rate values were determined as time averages of values provided by PIPE simulations. In order to test the validity of this method, we run atmospheric dispersion simulations with release rate steps (determined again as averages) of 1 s duration each, for some trials. Results of maximum DSL for the concentrations considered gave differences of