SciELO - Scientific Electronic Library Online

vol.36 número4Thermodynamic approach for optimal design of heat and power plants: Relationships between thermodynamic and economics solutions índice de autoresíndice de materiabúsqueda de artículos
Home Pagelista alfabética de revistas  

Servicios Personalizados




  • No hay articulos citadosCitado por SciELO

Links relacionados

  • No hay articulos similaresSimilares en SciELO


Latin American applied research

versión impresa ISSN 0327-0793

Lat. Am. appl. res. v.36 n.4 Bahía Blanca oct./dic. 2006


Atmospheric dispersion of emissions due to leakages in pressurized natural gas ducts

F. C. M. M. Soares, O. Q. F. Araújo and J. L. Medeiros1

Escola de Química - Universidade Federal do Rio de Janeiro, Centro de Tecnologia, Bloco E - Ilha do Fundão - Rio de Janeiro - RJ - Brazil - 21949-900.

Abstract — A simplified approach is presented to the transient atmospheric dispersion of accidental releases of natural gas, originated by leakages in pressurized ducts on sea level. In this scenario, shut-off valves are used for instantaneous shutdown of tube operation, isolating the harmful inventory, whose transient release is simulated providing estimations of gas instantaneous atmospheric concentration. The analysis further covers the transient behavior inside the ducts through a leakage model and the occurrence of multiple ruptures, synchronized or not, with known spatial distribution. The time-space dispersion model employed accounts for: (i) atmospheric conditions, (ii) wind speed, (iii) transient conditions of gas release into the atmosphere, and (iv) plume rise. The process of continuous release is approximated by a finite sequence of pulses, known as "puffs". Discrete puffs have volume dependent on release intensity, which depends on the actual time instant, as inventory decreases due to emission.

Keywords — Atmospheric Dispersion. Multi-Source Plumes. Plume Rise and Leakage.


Considering the flammability risk of materials in off-shore processes, the pressures and flow rates involved, the atmospheric releases of natural gas from leakages (or failures) in a platform and/or in its neighbor pipeline system configure events of high risk, which, in case of pipeline rupture, can cause significant fatalities and damages to the environment. The prediction of the ensuing release rate and its variation with time are the two critical pieces of information required in assessing and quantifying the consequences associated with such failures and may dictate the survival time of the mechanical integrity of the production platform (Mahgerefteh et al., 2003). An example of an offshore incident related to gas release is the propagation of fire from the Tartan riser failure, in 1988 on North Sea. This leakage set the destiny of the Piper Alpha platform, generating 167 deaths, total loss of the platform and hundreds of millions of dollars in losses (Cullen, 1990).

On the occurrence of an emergency, the main strategy of protection against a leakage of light hydrocarbons consists on isolation of the damaged duct and its dangerous inventory through shut-off valves, while the gas is liberated to the atmosphere. This solution is difficult to put into practice due to interconnecting ducts and risers, being the failure of those, in general, the main factor defining the risk level. The crew proximity to possibly dangerous situations and its confinement on vessels bring extra need to risk assessments, such as: (i) transient behavior on leakage flow rates; (ii) response time of detection and shutdown systems; (iii) time-space distribution of the leakage plume around the rupture point; (iv) effects of gas transport in the atmosphere.

The present work focuses on the occurrence of leakages on pressurized gas pipelines (risers and ducts to mainland) on production platform adjacencies. Specifically, the occurrence of ruptures on natural gas pipelines is investigated. Therefore, physical effects of a particular event are approached, such as discharge flow rate, alterations on ducts network flow and atmospheric dispersion of the leakage plume around the rupture point. The dispersion model employed is based on numerical resolution of the Tri-dimensional Diffusion Equation, under turbulence and on semi-infinite domains (Medeiros, 1993). The model accounts for: (i) ample variety of atmospheric conditions, changing the dispersion coefficients in function of the Pasquill stability classes recommended by Briggs (1973) and E.P.A. (1995); (ii) wind velocity, admitted constant, parallel to the sea level and corrected to the discharge height (E.P.A., 1995); (iii) transient conditions of gas release into the atmosphere, such as velocity, flow rate, pressure, temperature and density (Medeiros, 1993); (iv) plume rise as a function of buoyancy force, momentum, discharge height, wind velocity and atmospheric conditions (E.P.A., 1989). The analysis contemplates the impact of shut-off valves location, and spacing among them, as well as the gas release to the atmosphere.


Atmospheric dispersion of hydrocarbon clouds depends strongly on how the atmospheric emission interacts with the wind, according to stability conditions. Dispersion conditions may vary drastically on the same instant, on the same location, as function merely of a small height difference with respect to the observed emissions - Fig. 1 (Turner, 2002).

Fig. 1. Distinct Dispersion Conditions
(source: Turner, 2002)

Therefore, the dominant factor in a time-space dispersion model of the leakage plume is the determination of stability conditions related to the simulated atmospheric emission. The main atmospheric conditions are divided into six different classes of

Pasquill stability according to the following codification: A= extremely stable; B= moderately stable; C= lightly stable; D= neutral; E= lightly unstable; F= unstable.

By means of Table 1, defined by Turner (1970), it is possible to classify which stability category is more suitable to the studied situation.

Table 1. Pasquill Stability Classes (source: Turner, 1970)


A. Gas Liberation under Leakage Model

The proposed model considers that safety devices (e.g., shut-off valves for isolation of the affected line) are actuated instantaneously after the leakage occurrence. Consequently, the amount of natural gas effectively liberated into the atmosphere is restricted to the material contained in between two shut-off valves that isolate the rupture point. The emitted volume is defined by the duct diameter and the spacing between these two valves. Given the gas density at the initial leak instant, the total mass and moles of the liberated inventory are determined. The emission presents transient behavior since the duct internal conditions vary as the gas inventory is depleted in the confined section of the pipe: temperature, pressure and density decrease (Medeiros, 1993). Hence, the gas is considered ideal, expanding isentropically, in accordance with Eqs. (1), (2) and (3):


where γ is the ratio between the heat capacity at constant pressure and the heat capacity at constant volume of the natural gas; ν is the molar volume of the gas; Pe is the duct pressure on the discrete instant e; T is the temperature and ρ the density on the indicated time. The emission mass flow rate is determined according to two distinct flow conditions - chocked and non-choked - as defined by the ratio of critical pressures expressed in Eq. (4) (Ramsbill, 1986), where indexes 1 and 2 refer, respectively, to the internal and external pressures:


The discharge flow is classified as non-choked if or choked if . At the non-choked flow, the mass flow rate, Qe, and velocity, uho, of the leaking gas are defined by Eq. (5) and (6), where: Cd is the discharge coefficient that characterizes the orifice friction, considered as 1.00 at full bore rupture and 0.80 for the remaining situations; Ah is the cross-section area.


The gas velocity (uho) is calculated by Eq. (6) until it reaches a maximum critical value corresponding to the sonic flow rate, uho,s, determined from the sound velocity given by Eq. (7). At choked flow, the flow rate is independent from external conditions:


B. Atmospheric Dispersion Model

Figure 2 shows the scheme of the dispersion model. The leakage plume occurs around a virtual source (height H), exactly above the rupture point, considered the real source (hs). There is a physical limitation: the ocean surface generates a reflection region that demands correction through the insertion of a source of the same capacity as the real source (hs), denominated the reflection source (-hs). The height difference between the virtual source and the real source is called the plume height. Due to the transient conditions in which the gas is liberated to the atmosphere, and being the dispersion a three-dimensional phenomena, transient results of concentration are generated for a fixed height plane of interest (z). This value may be the virtual source height (H), where the flammable gas concentration is higher, or at the platform, at a potentially critical height of risk.

Fig. 2. Dispersion Model Scheme

Plume Rise. The elevation phenomena occurs as a function of wind velocity, buoyancy force, gas momentum flux at the leakage, atmospheric conditions and discharge height (E.P.A., 1989). In this work, the height of real source emission point, hs, is located on the ocean surface. However, to avoid physical limitations generated at the boundary layer (i.e., ocean surface with zero wind velocity) a value of 30cm was used. The wind velocity, ur, is measured at the reference height, hr, equals to ten meters, and is corrected through Eq. (9) to the real source height, hs. The exponent w is established as a function of the Pasquill atmospheric condition as shown on Table 2 (E.P.A., 1995).


Table 2. Dispersion Parameters and Wind Correction (E.P.A., 1995)

In these equations, Fb is the buoyancy force and Fm is the inertial term, related to the gas emission momentum:


The plume rise, Δh, is calculated for each emission source j through empirical equations recommended by the U. S. Environmental Protection Agency (1989):


where the gravity acceleration (g) is set as 9.81m/s2, dh is the bore diameter, uho is the emission velocity and the temperature (T) indexes 1 and 2 refer to the gas temperature and to the atmospheric temperature, respectively. The difference between these two temperatures is defined as Δ T.

On the simulation of unstable to neutral atmospheric conditions corresponding to Pasquill classes A, B, C and D, if the buoyancy term, Fb, is less than 55N, critical Δ T, defined as Δ Tc, is calculated by means of Eq. (12) where α is the rupture angle with respect to the horizontal plane. For Δ T > Δ Tc, the plume rise is calculated by Eq. (14) and, otherwise, by Eq. (15). At these same atmospheric conditions, for Fb higher than 55N, Δ Tc is calculated by Eq. (13) and the plume rise through Eq. (15) and (16), being the first appropriated to Δ T Δ Tc, and the second to Δ T > Δ Tc.

To simulate the dispersion under stable atmospheric conditions, E and F Pasquill classes, it is necessary to determine the square of the Brunt-Vaisala frequency, s, using Eq. (17). This quantity depends on the ambient temperature gradient as a function of height dθ /dz. This gradient is estimated as 0.02K/m for class E and 0.035K/m for class F (E.P.A., 1995). At these stable atmospheric conditions, Δ Tc is determined through Eq. (18). In case of Δ T > Δ Tc and α less or equal to 90 degrees, the plume rise depends on the corrected wind velocity, u, and the critical wind velocity, uc, defined by Eq. (19). If u uc, the plume rise is determined by Eq. (20), otherwise, through Eq. (21). On the other hand, in case of Δ T Δ Tc and/or α larger than 90 degrees, the plume rise is calculated using Eq. (22).


Discrete Clouds. The mass and volume of the leaking gas are quantified and discretized on time through clouds materialized instantaneously, on the exact instant of liberation, with the shape of a parallelepiped. After created, these clouds are immediately liberated to the wind convective drag. As they move in the wind direction, the clouds deform smooth and gradually by the action of mass transfer diffusion in the lateral (y), vertical (z) and axial (x) directions. The x axis is associated to the wind direction and orientation.

At each time interval e (te-1, te) previously established, a cloud (or puff) with index e is formed. This indication follows the temporal order of clouds formation, being e equal to 1 for the first cloud. Naturally, the larger the number of discrete clouds, the better will be the emissive process description of the source along time. The initial coordinates of the eth cloud emitted by jth source are defined by: (i) the atmospheric conditions, (ii) the quantity of gas released on the time interval of the cloud creation, and (iii) the spatial position of the emission source.

Considering, for each cloud, the Cartesian directions 1, 2 and 3 as referring, respectively, to its length (x), width (y) and height (z), the puff initial volume, Vp, is calculated through Eq. (23):


where a and b refer to the initial and final coordinates of an edge from the parallelepiped aligned to the indicated direction (1, 2 and 3). The puff coordinates are calculated by Eqs. (24) to (27):


where: rj is the Cartesian position vector of the virtual source j, varying as a function of the plume rise Δ hj; Δ t is the discrete time interval between cloud emissions of the emission source; Se is the puff total height, estimated with the mass leakage flow rate, Qe; u is the wind velocity; ρ 2 is the gas density at external conditions; Ryz is an approximation of the lateral to vertical size ratio, from Table 2, and defines the discrete cloud width as a function of the puff total height, according to the Pasquill classes (E.P.A., 1995).

The phenomena of discrete clouds dispersion from emission source j, is modeled with a function of the source generation (fj) defined by Eqs. (28) and (29), where δ is the unit impulse function, H(ri,ai,e,bi,e) is the unit pulse function at the Cartesian axis i on the [ai,e,bi,e] interval, and Cg is the gas concentration at the emission source, obtained from ideal gas equation at atmospheric conditions - Eq. (30). These relationships express that the continuous emission of source j is substituted by an intermittent operation where, at each instant te-1 (e=1,2,, a cloud with mean properties related to the time interval Δ t or e (te-1, te) is created and freed to disperse as it is dragged by the wind (Medeiros, 1993).


Transient Dispersion Model. Medeiros (1993) developed an approximated solution to the Tridimensional Equation of Diffusion with pollutant generation term and convective wind transport, v - Eq. (31) - in order to determine the concentration field C(x,y,z,t). In this expression, f(x,y,z,t) is the gas generation function (in volume), which relates to each source j (fj in Equation (28)) according to Eq. (32):


where nj is the total number of sources and N is the diffusive flux, calculated through the diagonal matrix of turbulent diffusivities, D, by analogy to the Fick's Law of molecular diffusion - Eq. (33). In Eq. (34), i=1,2 refer to the diffusivities in the lateral directions and i=3 relates to the vertical diffusivity. These diffusivities are calculated as a function of the position x at the wind convective direction and of the dispersion parameters, σ , found in Table 2 in accordance to the Pasquill stability classes (Briggs, 1973). k is a conversion factor used, equivalent to 1.0m/s. Further, the wind velocity vector v is considered constant, assuming only one component different from zero: v2 = v3 = 0; v1 = u.

To implement the resolution of problem defined by Eqs. (31) to (34), the following initial and boundary conditions were adopted: (i) absence of gas in the atmosphere at t = 0; and (ii) zero diffusive flux on the sea surface. The solution of this problem is given by Eq. (35), where Pi(ri,t,aij,bij,te-1,j) defines the contribution from cloud j to the distribution of gas concentration regarding the indicated direction i, in order to be: (i) null, in the case of t te-1,j; (ii) established, for t te-1,j, by means of Eq. (36). In this expression, Ai and Bi are matrices defined by Eq. (37) and (38).



The model was implemented in MATLAB and simulations were conducted to obtain information concerning: (i) the impact of shut-off valves spacing; (ii) the spatial configuration of emissions. The critical scenario defined by Soares et. al. (2004) corresponds to: (i) full bore rupture; (ii) strong wind straight towards the platform; (iii) Neutral - D - Pasquill stability condition. Results of simulations are presented for different risk situations due to three emission sources of natural gas. Ambient pressure and temperature were considered as 1.013bar and 25oC, respectively. Inside the duct, the pressure and temperature were initially at 20bar and 15oC. The duct diameter was 40cm. The height of simulation plane (height of interest) was the virtual source itself, where the highest flammable gas concentration occurs. Wind velocity of 6.0m/s at reference height was assumed.

The proposed spatial configuration was simulated for the neutral - D - category on three different cases that have, in common, the time interval between the failures at the three emission sources stipulated in four seconds. That is, at each four seconds, from the first leakage, considered instantaneous, occurs the next leakage, to the total of three leakages.

Figure 3 shows the three studied cases. In Case A, the distance between shut-off valves is defined as five kilometers. The ruptures, represented by circles, were positioned in series towards the platform, spaced ten meters. The platform was localized at 200m from the first leakage and, consequently, at 180m from the third source. The pipelines are schematically displayed as rectangles to the left of the platform equipments. For Case B, the ruptures were positioned in parallel, and the platform was 200m from the three sources. The distance between the valves was of ten kilometers and the emission sources were three meters apart. Case C refers to a situation in which the emission sources are placed in parallel and distant seven meters in the lateral direction (y axis). Each duct had an emission source with a shift of ten meters from the neighbor source in the axial (x) direction. The platform was located at 200m from the first source and the valves were distanced by five kilometers.

Fig. 3. Case Studies A, B and C

Results for Case A are shown in Fig. 4. The discharge ends before 63s have elapsed from the first leakage. Due to the serial positioning of the sources in Case A, the concentration that reaches the 200m platform front line at the virtual source height (z=84.6m) is very high (~35gmol/m3), generating the worst risk situation associated to the proposed scenario.

Fig. 4. Concentration Profile for Case A: z=84.6m; t=63s; Wind Velocity= 6m/s; Atmospheric Class D

Figure 5 is related to Case B after 63s: it can be noticed that leakage effect remains and, even though the emission sources were laterally distanced by three meters, there still is a non-zero gas concentration. At the 200m platform front line, the concentration is superior to 20gmol/m3. Additionally, the longer distance between the shut-off valves (10km) has practically doubled the leakage time to 120s, aggravating the associated risk.

Fig. 5. Concentration Profile for Case B: z=84.6m; t=63s; Wind Velocity = 6m/s; Atmospheric Class D

Case C (Fig. 6) may be considered the safest one in reason of the small distance between the valves (smaller gas inventory) and the seven meters of lateral distance, which avoids the concentration build up by additive emissions. The concentration that reaches the 200m at the virtual source height is less than 15gmol/m3 and the emission ends before 63s.

Fig. 6. Concentration Profile for Case C: z=84.6m; t=63s; Wind Velocity= 6m/s; Atmospheric Class D


This work presents a mathematical model for risk analysis of offshore production fields associated to the atmospheric dispersion of natural gas from leakages and ruptures at pipeline networks.

The model was implemented in MATLAB environment and describes the time-spatial distribution of gas concentration originated by one or more sources, in order to: (i) identify the worst scenarios; (ii) formulate solutions that effectively contribute to the minimization of risks. In addition, three source configurations of multiple leakages were simulated, which characterize a multi-source problem. The worst situation corresponded to ruptures occurring in serial arrangement.

1. Briggs, G. A. Diffusion Estimation for Small Emissions. NOAA ARL Rep. ATDL-106, 1973.         [ Links ]
2. Cullen, W. D. The public inquiry into the Piper Alpha Disaster. London Energy Dep., HMSO,1990.         [ Links ]
3. E.P.A. The Offshore and Coastal Dispersion Model - Volume I. Rep. No, A085-1, 1989.         [ Links ]
4. E.P.A. User's Guide for the ISC3 Dispersion Model - Volume II. EPA-454/B-95-003, 1995.         [ Links ]
5. Mahgerefteh, H., Oke, A., Economou, I., and Rykov, Y. A transient outflow model for pipeline puncture. Chemical Engineering Science 58, 4591-4604, 2003.         [ Links ]
6. Medeiros, J. L. Simulação Transiente de Dispersão Atmosférica de Poluentes Gasosos. Anais do XIV C. Ibero Latino-Americano de M. Comp. em Eng., IPT, S. P., p. 1120-1129, 1993.         [ Links ]
7. Ramsbill, P. K. Discharge Rate Calculation Methods. Saf. and Reliab. Dir., UKAEA., 1986.         [ Links ]
8. Soares, F. C. M. M., O. Q. F. Araújo and J. L. Medeiros, Simulação de cenários de risco associado à ruptura de dutos em plataforma Off-Shore. XV Cobeq, 2004.         [ Links ]
9. Turner, D. B. Workbook of Atmospheric Dispersion Estimates. EPA Office Air Prog. Rep., AP-26, 1970.         [ Links ]
10. Turner, D. B. Air Dispersion Modeling Basics. Carolinas Air Pollution Control Association, 2002 Fall Meeting Presentation, 2002.
        [ Links ]

Received: December 20, 2005.
Accepted for publication: July 15, 2006.
Recommended by Editor A. Bandoni.

Creative Commons License Todo el contenido de esta revista, excepto dónde está identificado, está bajo una Licencia Creative Commons