Stéphane Lafortune
email: stephane.lafortune@ineris.fr.

Zbigniew Pokryszka

Gaëtan Bentivegna

Régis Farret

From geochemical baseline studies to characterization and remediation of gas leaks: experiences and case studies of the French institute for risk management (INERIS)

Keywords: monitoring, gas, leak, surface monitoring, vadose zone, contingency plan, remediation plan

Involved in several research projects funded by France and European Union on carbon capture and storage, INERIS has worked for many years on the design, building, testing and field-deployment of gas monitoring devices in soil and subsurface on industrial sites. INERIS has developed a strong experience using self-designed devices dedicated to the characterization of CO 2 storages and analogues of natural or anthropogenic CO 2 -leaking systems. INERIS focuses on the monitoring of gas flux at surface and on the monitoring of gas migration in the subsurface using gas monitoring wells.

After many years of worldwide research done by many researchers and industrials, it's now time to adopt a "monitoring-ready" approach. That means to be ready to challenge industrial and stakeholders' expectations for efficient, field-deployable and cheap solutions to perform efficient gas monitoring in the soil and the subsurface.

We will here discuss the way INERIS performs surface and subsurface gas measurements (1) to achieve baseline studies on storage sites before injection starts, (2) to ensure monitoring in injection and post-injection phases and (3) to characterize leaking processes that impact surface and urban areas.

Introduction

Gas and water geochemistry has a major role to play in the surface and subsurface monitoring of storage sites. See for instance the concern that grew up from 2004 to 2011 about a potential CO 2 leak from the Weyburn-Midale Carbon Capture and Storage (CCS) project in Saskatchewan (Canada) that would be impacting a farm close to the oil-field and force the owners to move from their property [START_REF]The Kerr investigation: executive summary[END_REF].

Investigations that were performed to discuss the origin of the "gas anomaly" that was detected in Saskatchewan highlight the need of a good geochemical monitoring plan to help [START_REF]The Kerr investigation: executive summary[END_REF] to define baselines before any CO 2 injection starts, (2) to track gas leaks that could impact surface and (3) to re-establish public confidence in the feasibility of safe carbon storage.

Through our own experience and feedback, we will here discuss: The monitoring devices and the monitoring strategy that have to be considered to establish baselines before injection starts and to perform good surface and subsurface gas monitoring injection and postinjection phases; The characterization of CO 2 leaks and the understanding of the gas migration dynamics in the subsurface, introducing the buffer role that could play the unsaturated zone (or "vadose zone") as a very last barrier before CO 2 reaches surface and accumulates in buildings, tunnels, underground car parks, cellars… Discussions will be based on our involvement in research and industrial pilot sites but also on study of a leaking system in a former mining area.

Working on CO 2 leaking analogues (e.g. volcanic area) is very insightful and challenging in a scientific and technical point of view, because we face conditions that cannot be encountered on industrial sites unless in accidental context. But we must handle the analogy with care and be aware of its limits… by also working on industrial or research pilot sites.

Performing surface and subsurface gas monitoring above CO 2 storages

INERIS is involved in the monitoring of an industrial CCS pilot site for several years. We have taken part in the baselines studies and take now part to gas monitoring during injection phase, by: Performing gas flux measurements at surface (soil/atmosphere interface) [2;3]; Characterizing gas migration processes in the subsurface (~100 meters depth) through gas concentration measurements in the gas phase of a monitoring well [START_REF] Lafortune | First steps in coupling continuous carbon isotopic measurements with already proven subsurface gas monitoring methods above underground carbon dioxide storage sites[END_REF].

Surface monitoring using gas flux measurements

INERIS monitors CO 2 flux above the industrial storage site since September 2008 by using a patented accumulation chamber (fig. 1). Monitoring is performed quarterly on 35 locations. Numerous field measurement campaigns have been held today. Some of them have been held before injection starts (= baseline studies) [START_REF] Lescanne | The site monitoring of the Lacq industrial CCS reference project[END_REF].

As comparison, gas flux measurements have also been performed in other locations in France far from a CCS site. These measurements are done since several years on 12 different locations situated along 2 long geographic profiles (profile A: N-S and profile B: W-E, see map on fig. 2 right). They constitute a database giving a range for CO 2 natural fluxes emitted at surface by (micro)biological and pedogenetic activities in soil that are usually observed in France all over the country for different land use contexts (forest, meadow, agricultural field... see fig. 1 left). This range is from 0 to 18 cm 3 .m -2 .min -1 (0 to 54 10 -8 kg.m -2 .s -1 ). Results show a strong relation between CO 2 flux value and seasonal weather conditions. Atmospheric temperature seems to be a key factor controlling gas emission at soil-atmosphere interface (fig. 2).

Any other significant contribution to CO 2 emission at soil-atmosphere interface such as a CO 2 leak from a deep storage would increase the quantity of gas that is emitted by biological and pedogenetic activities in soil and would then be theoretically detected. It is a direct and very flexible monitoring method that can detect any significant anomaly in time as long as measurements are regularly performed and statistically workable; It is a light and fast monitoring method to quickly confirm or remove any doubt concerning the occurrence of a gas leak from storage at surface; It is a method that can be easily coupled with other measurements performed in soil such as gas soil measurements or isotopic analysis of the emitted gas; It is a helpful communication tool to establish public confidence in the safety of CCS project (this is a direct method which is visible and whose physics is easily understandable).

Subsurface monitoring using dedicated wells

Before reaching surface, gas leak from CO 2 storage will migrate through upper geological formations. Knowledge of underground residence time and processes that can impact this migration is fundamental to determine the potential consequence of a leak at surface.

INERIS works for several years on the characterization of gas migration processes in the unsaturated shallow formations (vadose zone) because these formations are the last barrier before gas reaches surface (see section 3 for a concrete example). As field case study, a monitoring well has been drilled close to the injection well of the CCS pilot site. This well was drilled down to about 100 meters depth and has reached water table at about 50 meters depth. Before reaching surface, a CO 2 leak from storage would migrate through the aquifer formation and eventually be dissolved in water. Degassing processes at water table could then occur and lead to CO 2 emissions in the vadose zone and at surface (fig. 3).

INERIS is monitoring gas concentrations and carbon isotopic signature in this well since several years. Atmospheric parameters (atmospheric temperature and barometric pressure) are also monitored. Measurements began before injection starts.

Aims are to better understand gas exchange dynamics between saturated and unsaturated zones and detect any potential anomalous gas migration towards surface.

Monitoring gas migration in subsurface has to be taken into account when drawing monitoring plan of CCS site because dynamics of underground gas exchanges can have a strong effect on gas emission at surface.

In some extent, underground porous media can constitute buffer volumes that accumulate gas (in dissolved or gaseous phases. When a change of the weather conditions occurs (for example a change of the atmospheric temperature but also of the barometric pressure as we will see in the section 3), the gas is released and can migrate towards surface.

Characterizing CO 2 emissions at surface in order to propose adapted remediation solutions

Since 2009, INERIS conducts research to understand production and migration processes of CO 2 in the vadose zone above a former underground mining area.

CO 2 emissions have been detected in urbanized areas at surface or in cellars. These emissions can expose people to O 2 -defficient and CO 2 -enriched air which can be hazardous (risk of unconsciousness or asphyxia) [START_REF] Lafortune | Understanding O2-deficient and CO2-enriched gas production and migration in the subsurface above a coal post-mining area through in situ gas monitoring and modelling[END_REF].

Even if there is no underground storage in this area, the context is very similar to what industrials could face in case of gas leak from CO 2 storage. We will see that feedback from this case study is very useful to help to draw remediation and contingency plans for CO 2 geological storages.

When gas emissions are detected at surface, a fast characterization has to be done to assess risks to people and if needed suggest emergency solutions.

This characterization should aim (1) to determine gas production mechanism or gas origin, (2) to identify migration pathways (diffuse or located emissions) and (2) to describe dynamics of the emissions (driving force and time variations). Results will influence content of remediation and contingency plans. Special attention should be paid to administrative authorities, inhabitants and press that can put pressure on people working on the characterization emissions to find out results and to design solutions. All measurements that have been performed by INERIS at surface and in the subsurface are listed and described in table 1.

Table 1. Measurements performed to understand production and migration processes of CO2 in the vadose zone of a former underground mining area

Aim

Measurement Protocol Determination of gas production mechanism or gas origin

CO2 concentration

Measuring CO2 and major gas concentrations in atmosphere of close spaces at surface (e.g. cellars) or in soil gas to detect any change of abundance.

Isotopic analysis

Measuring the carbon isotopic signature of CO2 to determine the gas origin [START_REF] Lafortune | First steps in coupling continuous carbon isotopic measurements with already proven subsurface gas monitoring methods above underground carbon dioxide storage sites[END_REF].

Identification of migration pathways

Gas flux at soil-atmosphere interface Performing CO2 flux measurements at surface using an accumulation chamber in order to map emissive areas.

Volume activity of radon

Monitoring volume activity of radon in gas emissions at surface to have clues on migration pathway and duration.

Description of dynamics of gas emissions

Gas pressure

Monitoring pressure gradient between underground and surface with pressure sensors to see if it can be a driving force of gas emissions.

Gas temperature

Monitoring temperature gradient between underground and surface with temperature sensors to see if it can be a driving force of gas emissions.

Gas flow

Monitor gas flow direction at highly emissive points (such as cracks, faults, wells…)

Results show that CO 2 can be produced by a geochemical process in the capillary fringe between saturated and unsaturated zones (at ~50 meters depth, see fig. 4). CO 2 accumulates in vadose zone and migrates later towards surface in low barometric conditions (sudden barometric pressure decreases are here the main driving force, see fig. 5). CO 2 preferentially migrates through the numerous cracks that are generally present on the whole extent of the former mining area. Usually hidden/covered by soil layer, cracks are not easily detectable unless gas anomalies have been already detected above them. Contingency plan drawn by local authorities implied design of efficient ventilation in close spaces (e.g. cellars) where CO 2 accumulations have been detected. Filling of some identified cracks has also been performed to stop gas emissions in the most sensitive areas (e.g. cracks very close to houses, emissions with high CO 2 enrichment or O 2 depletion…).

INERIS recommends for mid/long-term remediation plans to design buildings so as they are not permeable to gas influx and by isolating first floors from cellars.

Foundations of new buildings should be gas proof to stop any penetration of gas emitted by soil. Gas tightness of old buildings should be improved by filling any potential pathway for gas (as cables and pipes passage systems or cracks in walls or floor).

Gas influx in old buildings should also be controlled by creating a pressure gradient between buildings and ground.

Feedback from this case study illustrates the fact that the vadose zone can be seen as a buffer volume which may accumulate gas migrating underground. This gas can be flushed later when meteorological conditions favor gas emissions at surface.

Because the vadose zone is the very last natural barrier before a gas leak reaches surface, it must be correctly monitored to detect any potential accumulation and prevent hazard to people. The monitoring of the vadose zone has to be taken into account when designing monitoring plans of CO 2 storages.

Since mid-2012, INERIS takes part to a new research project called "CIPRES" (Study of potential impacts of CCS on water resource) funded by the French Research Agency (ANR). INERIS is presently building a pilot site to simulate a leak of CO 2 in subsurface (fig. 6). Dissolved CO 2 will be injected in a shallow aquifer at about 20 meters depth. CO 2 migration in the saturated and unsaturated zones will be monitored through observations wells. This project will aim to better characterize degassing process at water table and also dynamics of gas migration and accumulation in the vadose zone. All techniques described in table 1 have been used by INERIS several times in different contexts related to gas monitoring. All of these techniques have passed validation tests at lab, are field-deployable and ready to monitor right now CO 2 storage sites. They can be used to initially characterize storages sites during baseline studies and to monitor them during exploitation and post-exploitation phases. They can also be used to characterize anomalous gas emissions at surface (potential gas leak from the storage).

Conclusion

Working since 20 years on designing and testing monitoring techniques, INERIS has a range of tools ready to monitor right now CO 2 storages before injection starts (baseline studies) and during exploitation and post-exploitation phases.

In case of suspicion of a gas leak from CO 2 storage, these tools can also be helpful to draw contingency and remediation plans by characterizing:

Gas emissions at surface to assess risks to people by crossing dangerousness of emissions with safety stakes; Gas migration in the subsurface to define migration process by determining what are the driving force and migration duration.

After many years of worldwide research, it's now time to answer industrial and stakeholders' expectations for efficient, field-deployable and cheap solutions to perform monitoring of storage sites. Panel of monitoring tools proposed here by INERIS can help designing surface and subsurface monitoring plans for CO 2 storage projects that are presently initiated or will be launched in a close future.

Fig. 1 .Fig. 2 .

 12 Fig. 1. Performing gas flux measurements at soil-atmosphere interface using an accumulation chamber

Fig. 3 .

 3 Fig. 3. Performing gas monitoring in the subsurface using a monitoring well

Fig. 4 .Fig. 5 .

 45 Fig.4. Plot of carbon isotopic signature of CO2 gas versus inverse of CO2 gas concentration (monitored gas appears to be the result of a mixing between atmospheric gas and source gas enriched in CO2 and having a carbon isotopic signature which is characteristic of a carbonates dissolution process)

Fig. 6 .

 6 Fig. 6. Pilot site of CIPRES to study impacts of CO2 leak on subsurface water resource and CO2 migration and accumulation processes in vadose zone

Acknowledgements

Authors thank research projects CIPRES and SENTINELLE (supported by the French National Research Agency, ANR) and GEODERIS (the French public expert for risks related to post-mining) for funding this work.


R M Woolley
email: r.m.woolley@leeds.ac.uk

M Fairweather
email: m.fairweather@leeds.ac.uk

C J Wareing
email: c.j.wareing@leeds.ac.uk

S A E G Falle

C Proust

J Hebrard
email: jerome.hebrard@ineris.fr

D Jamois

Experimental Measurement and Reynolds-Averaged Navier-Stokes Modelling of the Near-Field Structure of Multi-phase CO 2 Jet Releases

Keywords: CCS, multi-phase flow, experimental measurement, mathematical modelling, accidental releases, atmospheric dispersion

The deployment of a complete carbon capture and storage chain requires a focus upon the hazards posed by the operation of CO 2 pipelines and intermediate storage vessels, and the consequences of accidental release. The aim of this work is the construction of a computational fluid dynamic model capable of accurately representing the complex physics observed in such a release, essential if dispersion phenomena are to be accurately predicted.

The interacting thermo-physical processes observed include those associated with the rapid expansion of a highly under-expanded jet, leading to an associated sonic flow structure. In such a release, it is also possible for three phases to be present due to the expansion and subsequent Joule-Thomson cooling, and a suitable equation of state is required to elucidate a system's composition. The primary objective of this work is the consideration of these physical processes, and their integration into a suitable numerical framework which can be used as a tool for quantifying associated hazards. This also incorporates the validation of such a model using data available in the literature and also using that recently obtained, and presented here for the first time. Overall, the model has provided an excellent level of agreement with experimental data in terms of fluid and sonic structure and temperature measurements, and good agreement with respect to composition data.

INTRODUCTION

Carbon capture and storage (CCS) refers to a set of technologies designed to reduce carbon dioxide emissions from large point-sources of production such as coal-fired power stations and other industrial facilities, in order to mitigate greenhouse gas release. CCS technology involves capturing CO 2 and then storing it in a suitable semi-permanent facility such as naturally formed saline aquifers or depleted oil wells, instead of allowing its release to the atmosphere where it contributes to climate change.

The technological requirements for the safe transportation of large quantities of liquids and gases at high pressures have been established over a number of years. Put into practice, the technology is proven for a number of hazardous substances including combustibles and toxins. Now however, with the advent of large-scale carbon capture and storage projects, attentions have turned to the safe transportation and handling of dense-phase carbon dioxide (CO2PipeHaz, 2009). The physics observed during the high-pressure release of combustibles such as natural gas are well studied and their behaviours relatively well understood. However, CO 2 poses a number of dangers upon release due to its more unusual physical properties, the physics of which have not yet been fully elucidated. Currently, the work being undertaken in the CO2PipeHaz project (CO2PipeHaz, 2009) is pivotal to quantifying all the hazard consequences associated with CO 2 pipeline failure, forming the basis for emergency response planning and determining minimum safe distances to populated areas. Such pipelines are considered to be the most likely method for transportation of captured CO 2 from power plants and other industries prior to subsequent storage, and their safe operation is of paramount importance as their inventory is likely to be several thousand tonnes.

CO 2 is a colourless and odourless asphyxiant which is directly toxic if inhaled in air at concentrations around 5%, and is likely to be fatal at concentrations around 10% [START_REF] Niosh | Documentation for Immediately Dangerous to Life or Health Concentrations (IDLHs) for Carbon Dioxide[END_REF]. Liquid CO 2 has a density much greater than water, but has a viscosity of magnitude more frequently associated with gases, and these properties make the transport of CO 2 an economically viable and attractive proposition. However, preliminary calculations and experimental evidence indicate that, due to it possessing a relatively high Joule-Thomson expansion coefficient, the rapid expansion of an accidental release may reach temperatures below -180 K. Due to this effect, solid formation following a pipeline puncture or rupture is to be expected, whether directly from liquid or via a vapour-solid phase transition with the passing of the system through the triple point (216.6 K at 5.11 atm). Additionally, CO 2 sublimes at ambient atmospheric conditions, which is behaviour not seen in most other solids. This is an important consideration when assessing the hazards posed by accidental releases, as CO 2 is denser than air, and an evaporating pool formed by liquid rain-out or a subliming solid deposit could cause a density-driven flow of high CO 2 concentration at ground level, in addition to the hazard posed by the associated gaseous release.

The developments presented in this paper describe a novel multi-phase discharge and dispersion model capable of predicting both the near-and far-field turbulent fluid dynamics and phase-transition phenomena associated with accidental CO 2 releases. As previously discussed, predicting the correct thermodynamic phase during the discharge process in the near-field is of particular importance given the very different hazard profiles of CO 2 in the gas and solid states. The modelling of CO 2 fluid dynamics therefore poses a unique set of problems, and the theoretical developments presented in this paper go some way to elucidating the observed physics. Previous works concerned with the near-field modelling of CO 2 releases are well covered in a recent review [START_REF] Dixon | Validation of Dispersion Models for High Pressure Carbon Dioxide Releases[END_REF], and should be consulted for other recent developments.

Also presented are the results of a series of experimental measurements of large-scale jetreleases of CO 2 , representative of pipeline punctures under various conditions. This series of experiments represents the most up-to-date data set available, describing the temperature and concentration fields of such under-expanded flows. Model validations have also been undertaken using the experimental data described, with shortcomings of the mathematical model elucidated through such comparisons, and suggestions for further developments presented.

It should also be noted that the modelling work presented here has further possible applications, and not solely within the field of CO 2 pipeline safety. Although the techniques used for the modelling of sonic jets are widely reported upon, the representation of the threephase CO 2 expansion and the associated thermo-physical phenomena is a novel approach. This could be applied in a number of areas of technology including those based upon the Rapid Expansion of Supercritical Solvents (RESS) processes including pharmaceutical, cosmetic, and speciality chemical industries. In these applications, the geometry of particles produced is determined by a number of factors including nozzle geometry, mass flow-rate, and pressure and temperature of CO 2 . Hence, an ability to model and predict fluid structures and particle distributions would be a great benefit in the design stage of such processes.

EXPERIMENTAL ARRANGEMENT

Figure 1 depicts a schematic of the rig used at INERIS for the experimental studies of largescale CO 2 releases. In the modelled region of the flow field, the instrumentation consists of fifteen radially distributed thermocouples, and three oxygen sensors distributed along the centre-line axis of the jet. The region used for the model validation extends to 5m from the release plane and is indicated by the shaded area of Figure 1. The choice of the modelling domain size is due to the model developments in this paper being concerned with the accurate representation of under-expanded, shock-laden, multi-phase jets, and the structure of their near-field. By 5m downstream of the release point, the jet has become self-similar in its properties, and has been at atmospheric pressure for a considerable distance. Hence, the modelling of this region does not require such specialist treatment.

Figure 2 is a photograph of the same rig, showing the external release point on the right of the picture, and the refilling sphere in-situ in front of the main vessel. The latter 2 cubic metre spherical pressure vessel is thermally insulated, and can contain up to 1000 kg of CO 2 at a maximum operating pressure and temperature of 200 bar and 473 K, respectively. It is equipped internally with 6 thermocouples and 2 high-precision pressure gauges as well as sapphire observation windows. It is connected to a discharge line of 50 mm inner diameter, with no internal restrictions. In total, the line is 9 metres long including a bend inside the vessel, plunging to the bottom in order to ensure that it is fully submersed in liquid CO 2 .

Three full-bore ball valves are installed in the pipe. Two are positioned close to the vessel and the third near to the orifice holder. The first valve closest to the sphere is a manual safety valve, and the two others are remotely actuated.

The vessel is supported by 4 load cells enabling a continuous measurement of the CO 2 content. The determination of the mass flow rate is performed within an accuracy of approximately 10%. The temperature is measured inside the vessel and immediately upstream of the orifice with 0.5 mm K-type thermocouples of accuracy better than 1 K. The static pressure is measured inside the vessel using a Kistler 0-200 bar instrument with an accuracy of ±0.1%, and immediately upstream from the orifice using a KULITE 0-350 bar instrument with an accuracy of ±0.5%. The vessel instrumentation is shown in Figure 3.

Various orifices can be and were used at the exit plane of the discharge pipe, and are all drilled into a large screwed flange. The thickness of this flange is typically 15 mm and the diameter of the orifice is constant over a length of 10 mm and then expanded with an angle of 45° towards the exterior. Figure 4 provides an example of such an orifice, whilst Figure 5 is a high-speed camera still of a typical release from a 9 mm nozzle. The discharge nozzle diameters used were 6, 9, 12 and 25mm in the six tests reported and studied here.

The field instrumentation consists of K-type thermocouples (0.5 mm diameter) arranged on vertical masts at varying distances from the orifice. The response time of these is approximately 1 second. Additionally, at each mast there is an O 2 sensor located on the centre-line of the jet.

A number of experimental configurations have been investigated during the project, but the most recent data are presented here. The details of six release scenarios, typical of CO 2 transport conditions, are presented in Table 1, which are used for the model validation.

MATHEMATICAL MODELLING

Governing Equations

Predictions were based on the solutions of the Reynolds-averaged, density-weighted forms of the transport equations for mass, momentum, two conserved scalars (CO 2 mass fraction and CO 2 dense phase fraction), and total energy per unit volume (internal energy plus kinetic energy), as described below by Equations 1, 2, 3, 4, and 5 respectively. This model is capable of representing a fluid flow-field comprising a mixture of CO 2 (vapour/liquid/solid) and air.

The equations employed in this study were cast in an axisymmetric geometry, but for brevity all subsequent equations are listed in their Cartesian tensor form.

( )

0 i i u t x ρ ρ ∂ ∂ + = ∂ ∂  (1) ( )  ( ) 0 i i j i j u j u u u p u u s t x ρ ρ ρ ∂ ∂ ′′ ′′ + + - -= ∂ ∂    (2)
( ) ( )

0 i t i i j u t x x x β ρβ ρβ µ   ∂ ∂ ∂ ∂ + - =   ∂ ∂ ∂ ∂         (3)
( ) ( )

0 i t i i j u s t x x x α α ρα ρα µ   ∂ ∂ ∂ ∂ + - -=   ∂ ∂ ∂ ∂        (4)  ( ) 0 i i i j t E i i j E S E p u u u u T s t x x x µ   ∂ ∂ ∂ ∂   ′′ ′′ + + - - -=       ∂ ∂ ∂ ∂       (5)
This equation set is closed via the prescription of the turbulence stress tensor (  i j u u ′′ ′′ ) as prescribed in Section 3.2.

Turbulence Modelling

Closure of the equation set discussed in Section 3.1 was achieved via the k ε turbulence model [START_REF] Jones | The Prediction of Laminarization with a Two-equation Model of Turbulence[END_REF], where  i j u u ′′ ′′ are the Reynolds stresses, and are modelled using this approach as:

 2 3 j k l i j ij t t k j l u u u u u k x x x ρ δ ρ µ µ   ∂   ∂ ∂ ′′ ′′ = + - +       ∂ ∂ ∂        (6)
Here, the turbulent or eddy viscosity is represented as a function of the turbulence kinetic energy by: 2 with 0.09

t k C C µ µ µ ρ ε = = (7)
which leaves the requirement for the solution of transport equations for the turbulence kinetic energy and its dissipation rate. Subsequently, these are modelled as:

( ) ( )

0 i t k i i i k k ku s t x x x ρ ρ µ   ∂ ∂ ∂ ∂ + - -=   ∂ ∂ ∂ ∂    (8)
and ( ) ( )

0 i i i i u s t x x x ε ε ε ρε ρε µ   ∂ ∂ ∂ ∂ + - -=   ∂ ∂ ∂ ∂    with 1.3 t ε µ µ = (9)
where the source terms k s and s ε are defined as:

k t s P ρε = - and
( )

1 2 t s C P C k ε ε ρε = - ( 10
)
with 1 2 1.4 and 1.94

C C = = (11)
and

     2 3 i i j i j t t t j j i i j u u u u u P k x x x x x µ ρ µ       ∂ ∂ ∂ ∂ ∂ = + - +           ∂ ∂ ∂ ∂ ∂         (12)
Although the standard k-ε model has been extensively used for the prediction of incompressible flows, its performance is well known to be poor in the prediction of their compressible counterparts. The model consistently over-predicts turbulence levels and hence mixing due to compressible flows displaying an enhancement of turbulence dissipation. A number of modifications to the k-ε model have been proposed by various authors, which include corrections to the constants in the turbulence energy dissipation rate equation [START_REF] Baz | Modelling Compressibility Effects on Free Turbulent Shear Flows[END_REF][START_REF] Chen | Computation of Turbulent Flows Using an extended k-epsilon Turbulence Closure Model[END_REF], and to the dissipation rate itself [START_REF] Sarkar | The analysis and modelling of dilatational terms in compressible turbulence[END_REF][START_REF] Zeman | Dilational Dissipation: The Concept and Application in Modeling Compressible Mixing Layers[END_REF]. Previous works by the present author (Fairweather andRanson, 2003, 2006) have indicated that for flows typical of those being studied here, the model proposed by Sarkar et al. [START_REF] Sarkar | The analysis and modelling of dilatational terms in compressible turbulence[END_REF] provides the most reliable predictions. This model specifies the total dissipation as a function of a turbulent Mach number and was derived from the analysis of a direct numerical simulation of the exact equations for the transport of the Reynolds stresses in compressible flows. Observations made of shock-containing flows indicated that the important sink terms in the turbulence kinetic energy budget generated by the shocks were a compressible turbulence dissipation rate, and to a lesser degree, the pressure-dilatation term.

In isotropic turbulent flow, the pressure-dilatation term was found to be negligibly small, and so it was proposed that the compressible dissipation rate take the form:

2 c t CM ε ε = (13)
where the constant C is set to unity to allow for the neglected pressure-dilatation term and t M is the turbulent Mach number. The application to the k-ε model is then made by modification to the source term of the turbulence energy evolution equation and to the turbulence viscosity as defined by Equations ( 14) and ( 15) respectively:

2 k t s M ρ ε = - (14) 2 2
(1 )

t t k C M µ µ ρ ε = + (15)
The turbulent Mach number is defined as:

( )

1 2 2 t k M c = ( 16
)
where c is the local speed of sound.

Figure 6 demonstrates the effect of these modifications upon the axial centre-line velocity profile predictions of a highly under-expanded air jet, reported by [START_REF] Donaldson | A study of free jet impingement. Part 1. Mean properties of free and impinging jets[END_REF]. The standard k-ε model is clearly too dissipative, leading to an early decay of the compression/decompression cycle. The compressibility corrected model, although remaining unable of describing the magnitude of the velocity at the peak of the second decompression cycle, evidently improves agreement with experimental observation by reducing the spreading rate of the jet.

Fluid Properties and Species Transport

Non-ideal Equation of State

The Peng-Robinson equation of state [START_REF] Peng | A New Two-Constant Equation of State[END_REF]) is satisfactory for predicting the gas phase properties of CO 2 , but when compared to that of [START_REF] Span | A New Equation of State for Carbon Dioxide Covering the Fluid Region from the Triple-Point Temperature to 1100 K at Pressures up to 800 MPa[END_REF], it is not so for the condensed phase. Furthermore, it is not accurate for gas pressures below the triple point and, in common with any single equation, it does not account for the discontinuity in properties at the triple point. In particular, there is no latent heat of fusion.

Span and Wagner [START_REF] Span | A New Equation of State for Carbon Dioxide Covering the Fluid Region from the Triple-Point Temperature to 1100 K at Pressures up to 800 MPa[END_REF] give a formula for the Helmholtz free energy that is valid for both the gas and liquid phases above the triple point, but it does not take account of experimental data below the triple point, nor does it give the properties of the solid. In addition, the formula is too complicated to be used efficiently in a computational fluid dynamics code. A composite equation of state has therefore been constructed to determine the phase equilibrium and transport properties for CO 2 . The inviscid version of this model is presented in detail elsewhere [START_REF] Wareing | A Composite Equation of State for the Modelling of Sonic Carbon Dioxide Jets[END_REF] and the method reviewed here is now extended for the turbulent closure of the fluid equations detailed in the previous section. In this, the gas phase is computed from the Peng-Robinson equation of state [START_REF] Peng | A New Two-Constant Equation of State[END_REF], and the liquid phase and saturation pressure are calculated from tabulated data generated with the Span and Wagner [START_REF] Span | A New Equation of State for Carbon Dioxide Covering the Fluid Region from the Triple-Point Temperature to 1100 K at Pressures up to 800 MPa[END_REF]

( ) ( ) ( ) ( ) a T RT P v b v v b b v b δ = - - + + - ( 17
)
where T is the temperature, R the Universal gas constant, and v is the molar volume. The parameters a and b are given to be:

( ) 2 0.45724 and 0.07780

crit crit crit crit RT RT a b p p = = (18) with ( ) 2 2 1 0.37464 1.54226 0.26992 1 crit T T δ ω ω     = + + - -             (19)
In Equations ( 18) and ( 19), the critical temperature, pressure, and acentric factor associated with CO 2 are given as: 6 7.3773 10 , 304.1282 and 0.228

crit crit p Pa T K ω = × = = ( 20
)
The internal energy in the gas phase is computed using the following expression for the specific heat at constant volume: The internal energy of the liquid is taken from Span and Wagner [START_REF] Span | A New Equation of State for Carbon Dioxide Covering the Fluid Region from the Triple-Point Temperature to 1100 K at Pressures up to 800 MPa[END_REF], except that 4 1

1.4422 10 J kg - × is added to the values in order to ensure that differences between the gas and liquid internal energies on the saturation line in the model are in agreement. The solid internal energy is described as: 2.04932 10 J kg - ×

( ) ( ) ( ) 2 2 3 3 5 1 4.04533 10 2 3 t t t T T T T U A T T B C J kg - - - = - × + -+ + (
. Note that the pressure dependence and the difference between the internal energy and the enthalpy for the solid are ignored since these are negligible.

To calculate the solid density, the same approach as Witlox et al. [START_REF] Witlox | Modelling of discharge and atmospheric dispersion for carbon dioxide releases[END_REF] is used, and expressed as: again based on property information from the DIPPRR 801 Database. From Liu [START_REF] Liu | Compression and phase behavior of solid CO 2 to half a megabar[END_REF], the solid sound speed at atmospheric pressure and 296.35 K is 1 1600 m s -and it is assumed that this is independent of temperature and pressure. Note that the results given below are extremely insensitive to the solid density and sound speed.

The saturation pressure above the triple point is taken from [START_REF] Span | A New Equation of State for Carbon Dioxide Covering the Fluid Region from the Triple-Point Temperature to 1100 K at Pressures up to 800 MPa[END_REF]. Below the triple point, they give the following empirical formula:

( )

1.9 2.9 Figure 7 shows the internal energy of the gas and condensed phases on the saturation line.

1 2 3 exp 1.0 1.0 1.0 t s t t t t T T T T p T p a a a Pa T T T T               = - + - + -                       (
The transition from liquid to solid has been smoothed over 4 K with a hyperbolic tangent function centred on the triple point. This has been done for computational reasons in order to ensure the function and its differentials are smooth.

When working with an equation of state, it is convenient to use the Helmholtz free energy (H) in terms of temperature and molar volume as all other thermodynamic properties can be readily obtained from it. For an ideal gas, the Hemholtz free energy per mole is given by: ( )

0 0 1 ln ln 1 i T v H H RT T v γ       = = - +       -       (26)
where 0 T and 0 v are reference states at which the ideal equation of state is valid. In the present non-ideal case, we can use the standard relation describing pressure:

T H p v ∂   = -  ∂   (27)
to obtain the Hemholtz free energy from any equation of state as:

( ) ( ) 0 0 , , v i v H T v H T v p d v = -∫ (28)
With a definition for the Helmholtz free energy, the entropy is obtained from the relation:

v H S T ∂   = -  ∂   (29)
and the internal energy is then given by:

( )

, H U T v H TS H T T ∂   = + = -  ∂   (30)
The local speed of sound, required by Equation ( 16) to construct the turbulent Mach number can now be derived. Using standard relations, it can be shown that:

1 S T v T v p p p S S v v T v T - ∂ ∂ ∂ ∂ ∂           = -           ∂ ∂ ∂ ∂ ∂           (31)
and:

1 1 2 2 2 S S p p c v v ρ       ∂ ∂   = -         ∂ ∂         (32)

Homogeneous Equilibrium Model

In a homogeneous equilibrium model (HEM), all phases are assumed to be in dynamic and thermodynamic equilibrium. Id est they all move at the same velocity and have the same temperature. In addition, the pressure of the CO 2 vapour is assumed to be equal to the saturation pressure whenever the condensed phase is present. The pressure of the condensed phase CO 2 is assumed to be equal to the combined pressure of CO 2 vapour and air (the total pressure). These assumptions are reasonable provided the CO 2 liquid drops or solid particles are sufficiently small so that they do not interact with the gas-phase turbulence

As the code employed in this work is conservative, it works with the total energy per unit volume as given by Equation 5, and also computes the total mass fraction of CO 2 ( β ) via Equation 3. In order to integrate the conservation equations, it is also necessary to calculate the total pressure, temperature, total mass fraction of CO 2 , density of CO 2 vapour, density of air, and density of condensed phase CO 2 , from the total density, U , and β .

In unit volume, the mass of the condensed phase CO 2 is then:

l m αβρ = (33)
which means the volume of the condensed phase CO 2 is:

l l V αβρ ρ = ( 34
)
The volume of the combined vapour and air is then:

1 1 g l l V V αβρ ρ = -= - (35)
Since the masses of the vapour and air per unit volume are:

( ) ( ) 1 and 1 v a m m β α ρ β ρ = - = - (36)
their densities are then:

( )

1 1 v v g l m V β α ρ ρ αβρ ρ - = =   -     (37) and
( )

1 1 a a g l m V β ρ ρ αβρ ρ - = =   -     (38)
Since the CO 2 vapour is in equilibrium with the solid/liquid CO 2 , the following holds:

( )

, v s v m p T p T ρ   =     (39)
where ( )

, v p T
ρ is the pressure given by the equation of state, Equation (17). In regions where there is significant mixing, one can use the ideal equation of state for the CO 2 vapour and:

( )

, v v v R T p T m ρ ρ = (40)
The total pressure is then given by:

a v a v p RT m m ρ ρ   = +     (41)
and the total internal energy by:

( ) ( ) ( ) ( ) ( ) 1 1 , 1 1 l v v a a U RT U T m m β α β αβ ρ γ γ   - - = + +   - -   (42)
where ( )

, l U T
ρ is the internal energy per unit mass of condensed phase CO 2 . The solid density is then determined from:

( )

, l l p T ρ ρ = (43)
which is obtained from the equation of state. Equations ( 33) to ( 38) are solved for T , p , and α using a Newton-Raphson iteration.

Homogeneous Relaxation Model

As previously discussed, the assumptions associated with the HEM are reasonable provided the CO 2 liquid drops or solid particles are sufficiently small. There are some indications that this will not be true, in particular for test calculations in which the release is from a nozzle with a diameter of the order of centimetres. Hence, the model was further developed as a homogeneous relaxation model, in that a relaxation time was introduced with respect to the transport of the dense phase. It was found that this extension to the model had significant effect upon predictions, and was required to effect the more accurate results presented in Section 4.

A full model requires the inclusion of discrete drops or particles, but it is possible to derive a simple sub-model for the relaxation to equilibrium in which the temperature relaxation is ignored and it is simply assumed that the condensed phase mass fraction is given by Equation This is consistent with the form of the evaporation/condensation rate given in Jacobson [START_REF] Jacobson | Fundamentals of Atmospheric Modelling[END_REF]. The relaxation time chosen to represent behaviour in the near-field of releases such as those considered herein was in the order of 10 -3 s and obtained by the assessment of the rate that the calculated CO 2 saturation pressure relaxed to the local vapour pressure. In post-shock regions of the flow, a relaxation time of the order 2.5 s was chosen, representing the non-equilibrium state of the condensed phase.

Discretisation and Mesh Adaption

Solutions of the time-dependent, axisymmetric forms of the descriptive equations were obtained using a modified version of a general-purpose fluid dynamics code referred to as MG, and provided by Mantis Numerics Ltd. Within this code, integration of the equations employed a second-order accurate, upwind, finite-volume scheme in which the transport equations were discretised following a conservative control-volume approach, with values of the dependent variables being stored at the computational cell centres. Approximation of the diffusion and source terms was undertaken using central differencing, and a Harten, Lax, van Leer (HLL) [START_REF] Harten | On Upstream Differencing and Godunov-type Schemes for Hyperbolic Conservation Laws[END_REF] second-order accurate variant of Godunov's method applied with respect to the convective and pressure fluxes. The fully-explicit, time-accurate method was a predictor-corrector procedure, where the predictor stage is spatially first-order, and used to provide an intermediate solution at the half-time between time-steps. This is then subsequently used at the corrector stage for the calculation of the second-order fluxes. A further explanation of this algorithm can be found elsewhere [START_REF] Falle | Self-similar Jets[END_REF].

The calculations also employed an adaptive finite-volume grid algorithm [START_REF] Falle | Body capturing using adaptive cartesian grids, Numerical methods for fluid dynamics[END_REF] which uses a two-dimensional rectangular mesh with grid adaption achieved by the successive overlaying of refined layers of computational mesh. Figure 8 demonstrates this technique. Where there are steep gradients of variable magnitudes such as at flow boundaries or discontinuities such as the Mach disc, the mesh is more refined than in areas such as the free stream of the surrounding fluid. Each layer is generated from its predecessor by doubling the number of computational cells in each spatial direction. This technique enables the generation of fine grids in regions of high spatial and temporal variation and, conversely, relatively coarse grids where the flow-field is numerically smooth. Hence, the code uses a hierarchy of grids, 0

N G G ⋅⋅⋅ such that if the mesh spacing is n x ∆ on grid n G then it is 2 n x ∆ on 1 n G + . Grids 0
G and 1 G cover the entire domain, but finer grids only exist where they are required for accuracy. Unlike codes such as FLASH [START_REF] Chicago | Flash Center for Computational Science[END_REF] that use patches of fine grid, refinement is on a cell-by-cell basis. The solution is computed on all grids and refinement of a cell on n G to 1 n G + occurs whenever the difference between the solutions on 1 n G -and n G exceeds a given error. Again, unlike FLASH [START_REF] Chicago | Flash Center for Computational Science[END_REF], each grid is integrated at its own time-step.

With respect to the establishment of grid-independent solutions, all calculations reported upon here were undertaken using refined grids in which the magnitude of fractional difference between the solution on the uppermost grid and a next level of possible refinement was below 0.5%. It has been ascertained through the study of numerous calculations, that this level of refinement provides solutions which purport negligible differences to those obtained on higher-resolution grids. In the case of the under-expanded air-jet calculations reported upon in Section 3.2, this resulted in fully-refined axi-symmetric meshes containing in the region of 130,000 nodes at 5 levels of refinement. Comparing this to previous similar studies, a non-uniform axisymmetric grid was reported to provide grid-independent solutions at a resolution of 126,000 nodes [START_REF] Birkby | Numerical predictions of turbulent underexpanded sonic jets using a pressure-based methodology[END_REF]; an adaptive technique, again used on an axisymmetric grid, reported similar success using 30,000 nodes [START_REF] Bartosiewicz | Numerical Investigations on Dynamics and Heat Transfer in a Turbulent Underexpanded Jet[END_REF]; and a full 3-dimensional model of a supersonic jet in a cross-flow was reported to be grid-independent using 1,544,098 cells [START_REF] Viti | Comparison of First and Second Order Turbulence Models for a Jet/3D Ramp Combination in Supersonic Flow[END_REF]. Hence, the required level of refinement used in the present work appears to agree with that noted by previous authors.

Additionally, it was found that the computational grids applied to the modelling of the CO 2 jets, required upwards of 500,000 nodes to achieve an independent solution. This was mainly due to the larger computational region which was required, extending to 5 metres downstream.

RESULTS AND DISCUSSION

The calculations required initial conditions of temperature, pressure, density, velocity, and dense-phase CO 2 mass fraction, which were obtained from isentropic decompression calculations of CO 2 . Initial conditions for the decompression were prescribed from experimental observation of conditions within the storage sphere, and at the head of the release pipe, which was assumed smooth. A sensitivity analysis was undertaken as to the effect of nominal pipe-wall roughness, and negligible effect was observed upon the results. It was however observed that the accuracy of the predictions was notably affected by the initially prescribed flow rate.

Figure 9 shows predictions of temperature profiles plotted against experimental data for tests 2, 3, and 4, at axial locations of 1, 2, 3, 4, and 5 m, observed on the centre-line of the jets.

The absence of data at the 1m position in test 2 is due to the freezing of the thermocouple and the consequent production of an erroneous signal. Additionally, Figure 9 compares predictions of O 2 molar concentration with experimental observation where data are available, namely centre-line locations in jets 2 and 3. These three tests were undertaken using the 6 mm nozzle, and data was only collected along the centre-line of the releases. The experimental sensor matrix was expanded in further tests to incorporate radial profiles, which will be subsequently discussed.

Predictions of axial temperatures of all three tests are in good agreement with observation, although an over-prediction of temperature in the near-field is observed in each. For tests 2 and 3, however, this over-prediction mostly lies within the limits of experimental error, but test 3 results show a slightly larger discrepancy at the one and two metre locations. This does however have the effect of bringing the far-field predictions at 2, 4, and 5 metres more into line with data. A slight change in gradient can be seen in the curve of predicted temperatures just before the axial location of 1 metre in all cases. This phenomena is to be expected, and is due to the system passing through the triple point as the temperature increases. The change in equilibrium from solid-gas to liquid-gas causes a change in the internal energy of the dense phase and hence effects a change in the local temperature. Additionally, predictions of O 2 concentrations along the centre-line of tests 2 and 3 are in good agreement with experiment.

An over-prediction of mixing within the first two metres of the release is translated into a slight over-prediction of O 2 fractions, but the models' performance remains meritorious. O 2 data for test 4 were not available. The over-dissipative solution is indicative of an incorrectly predicted mixing rate, and the k-ε turbulence model is well known to underperform in such a manner in compressible jets such as these. Although corrected according to the model of [START_REF] Sarkar | The analysis and modelling of dilatational terms in compressible turbulence[END_REF] there is the possibility of an anisotropic element of the Reynolds-stress tensor not being accounted for, and hence a second-moment turbulence closure is currently being incorporated within the model framework to address this issue.

Figure 10 shows predictions of radial temperature profiles plotted against experimental data measured along a vertical line through the release for tests 6, 7, and 8, at axial locations of 1, 2 and 5 m. The model qualitatively and quantitatively captures the thermodynamic structure of the sonic releases, and although there is a small discrepancy between the observed and predicted spreading rates in the very near-field, calculations lie within the accepted error range of the experimental data. Within the core of the jets, temperatures are seen to be slightly under-predicted when compared to experiment, except in the predicted inviscid region still present at 1 metre in tests 6 and 7, and 2 metres in test 8. It is possible that dense phase CO 2 is removed from the system due to such phenomena as agglomeration, which would affect the higher temperatures observed. Hence, recent developments of the model include the incorporation of sub-models for the distribution of solid and liquid particles within the flow, and it is expected that the effects of phenomena such as particle coagulation will have an impact upon the predicted temperatures. Also, the system may not be in equilibrium due to this, or the generated turbulence, which may cause the observed discrepancies. The peaks locate the centre of the jets, where pure CO 2 is present in both solid, and vapour phases. The solid and vapour are in equilibrium here, and hence the CO 2 vapour pressure is equal to the CO 2 saturation pressure. The local temperature is thus defined by this condition in the model. The total pressure in this region is near to atmospheric and hence a high concentration of CO 2 will affect a higher local temperature. Moving away from the centre-line, and into a region of mixing, the CO 2 fraction is reducing, leading to a reduction in vapour pressure and hence a lowering of the temperature. Once no solids remain, there is an inflection in the curve and the temperature is seen to increase with the mixing of the CO 2 vapour and ambient air. By 2m in tests 6 and 7 and by 5 metres in test 8, there is no dense phase predicted as remaining and the temperature profile is determined by the mixing of CO 2 with air at atmospheric pressure. It should be noted that this centre-line peak in predicted temperatures appears to be a feature of the homogenous model applied in these calculations, and this numerical phenomena is not reproduced in experimental observation, although further more detailed measurements would be useful in this regard. It is anticipated that the removal of the effect of the HEM assumption, and hence the formation of the troughs in the off-centre thermal profile, will go some way to bring predictions more in line with data. It is therefore expected that this will not be observed in calculations incorporating recently developed Lagrangian particle tracking models, which is the subject of further work being currently undertaken. 1metre in these investigations is however notably accentuated when compared to the previous experiments, and requires further investigation. According to [START_REF] Birch | The Structure and Concentration Decay of High Pressure Jets of Natural Gas[END_REF], the axial concentration decay of an under-expanded jet release can be prescribed as:

1 2 ps a g kd z n ρ η ρ   =     +   (44)
where the pseudo-diameter is given by: 0 0.587

ps d a p d d C p   =     (45)
where 0.587 is a CO 2 dependent constant. The axial decay constant k is given a universal value of 4.90, and d C a value of 0.85. This approach to the approximation of the mixing in tests 6, 7, and 8 was applied, and associated O 2 concentrations calculated. The crossed-circle symbols of Figure 11 show these values plotted against experimental observation and predicted results. Interestingly, the conformity of experimental observation and the prescription of [START_REF] Birch | The Structure and Concentration Decay of High Pressure Jets of Natural Gas[END_REF], decreases with increasing nozzle size and initial reservoir pressure. It is also notable that the experimental data are in close agreement with the theory of (Birch et al., 1984Birch et al., 1984) in the near-field, and conversely closer in agreement with predictions farther downstream. This is most notable in test 8 in which discrepancies at all three measurement points are large. Although not conclusive, these observations are indicative of the previously discussed possible shortcomings of the model.

CONCLUSIONS

A turbulent computational fluid dynamic model capable of predicting the near-field structure of high pressure releases of multi-phase carbon dioxide representative of those arising from an accidental pipeline puncture or rupture has been presented. This model incorporates a novel approach for the evaluation of the state variables of CO 2 , as recently presented by Wareing et al. (2012). Validation of this model has been undertaken against new experimental data also reported herein, and gathered as part of the CO2PipeHaz project.

Overall, the model developed has yielded an excellent level of agreement with measurements of the temperature characteristics of these jets, and good agreement with composition data in a number of comparisons.

It has been identified that two areas of improvement are required to ensure accurate representation of the complex physics observed in these release scenarios and developmental work is currently ongoing to address these issues. Firstly, a second-moment turbulence closure is required to represent the turbulence anisotropy, which is expected to correct errors in the predicted rates of mixing due to the two-equation model implemented. Secondly, the inclusion of particles within a Lagrangian framework is required to more accurately represent the thermophysical interactions between the phases. A model describing the motion of fluid and solid particles has been written into the code, and is currently under validation. Both droplet and particle populations are predicted through solutions of equations for the particle distribution function used for small particles, with a Lagrangian particle tracking routine employed for larger particles.

Additionally, it is clear from the predictions of these dense phase releases that significant solids are generated within the near-field of these jets, despite the release itself containing no dry ice. This is an important conclusion with respect to the future design of CO 2 pipelines and the consideration of the related hazards.

The safety issues surrounding the transport of CO 2 in pipelines are focussed on the questions of the toxicity of an accidental release and the fast decompression of a pipeline following a catastrophic extension of a crack initiated at a local puncture. Here, the case of pure CO 2 has been investigated, but the impact of the presence of annex gases and impurities was considered out of the scope of the present work. It is to be expected that a number of different impurities may be present in a captured CO 2 stream dependent upon the industrial source, and the presence of these additional compounds, some of which are more toxic than CO 2 (e.g. SO 2 , H 2 S, etc.), can drastically affect the temperature and pressure fields in the immediate vicinity of a local leak, eventually leading to fast crack propagation and total pipe failure.

Also, and depending on the interactions between the bulk CO 2 and impurities in the pipe, the dispersing cloud conditions could also be so strongly affected that its composition may differ significantly from the composition of the transported fluid, potentially leading to clouds that are more hazardous because of higher levels of toxics due to the mechanism of phase partitioning occurring just outside the pipe, very close to the leakage point.

The presence of impurities therefore has a notable effect on the phase behaviour of CO 2 , and hence alters the near-field characteristics of any accidental release. The introduction of small amounts of N 2 into high pressure CO 2 , for example, shifts the triple point towards higher pressures and lower temperatures. These variations modify the near-field of the dispersing jet in terms of its shock structure, temperature and pressure which in turn impact on the phase distribution of the dispersing CO 2 and impurities. It is therefore concluded that the next step in code development will be the extension of the methodologies discussed here to model the near-field structure and dispersion characteristics of releases of CO 2 mixtures containing impurities and trace elements typical of those found in CCS streams. This will include both a continuous phase fluid model, and as previously discussed, a discrete-phase model capable of representing the distribution and properties of solid and liquid particles in such releases. The latter is an important consideration given the unusual phase behaviour of CO 2 , and its sublimation from solid form at atmospheric conditions.

NOMENCLATURE

FIGURE CAPTIONS

TABLE CAPTIONS

Table 1. Parameters of the experimental releases.

.

  This expression has been obtained by fitting to the internal energy data in the DIPPR database (DIPPR, 2013) over the temperature range energy of the gas at the reference state used by Span and Wagner (STP) is then:

  been obtained from the DIPPRR Database tables along with the latent heat of fusion at the triple point of 5 1

  vapour pressure, s p the saturation pressure, and τ the relaxation time.

Figure 11

 11 Figure 11 depicts axial profiles of temperature predictions plotted against experimental data along the centre-line in the CO 2 releases of tests 6, 7, and 8. As previously discussed, the level of agreement between calculation and experiment is comparable for these three investigations. Also reflected is the observed centre-line under-prediction of temperature common to tests 6 and 7. Figure 11 also displays predictions of O 2 molar fraction plotted against experimental data on the centre-line of the same tests. Contrary to calculations of tests 2, 3, and 4, an under-prediction of temperature in the near-field now leads to an associated under-prediction of O 2 fraction at the same locations. This under-prediction at the distance of

Figure 1 .

 1 Figure 1. Schematic diagram of INERIS CO 2 release test rig including sensor configuration.

Figure 2 .

 2 Figure 2.General view of experimental rig, including filling sphere, main vessel and discharge pipe.

Figure 3 .

 3 Figure 3. Pressure vessel instrumentation.

Figure 4 .

 4 Figure 4. Example of orifice flange.

Figure 5 .

 5 Figure 5.High-speed camera still of a 9mm release.

Figure 6 .

 6 Figure 6.Predicted normalised axial centre-line velocity of a highly under-expanded air jet obtained using a standard and a compressibility-corrected k-ε turbulence model, plotted against experimental data.

Figure 7 .

 7 Figure 7. CO 2 internal energy predictions on the saturation line using the composite equation of state, showing gaseous and dense phases.

Figure 8 .

 8 Figure8.Adaptive mesh refinement grid mapped onto velocity predictions in the region of a Mach disc.

Figure 9 .

 9 Figure 9. Predictions (lines) of axial temperature and O 2 mole fraction profiles (not test 4) plotted against experimental data (symbols) for tests 2, 3, and 4.

Figure 10 .

 10 Figure 10.Predictions (lines) of vertical temperature profiles plotted against experimental data (symbols) for tests 6, 7, and 8, at axial locations of 1, 2 and 5 m.

Figure 11 .

 11 Figure 11. Predictions (lines) of axial temperature and O 2 mole fraction profiles plotted against experimental data (symbols) for tests 6, 7, and 8.

Figure 03 Figure 04

 0304 Figure 03

  Figure 05

  Figure 08

ACKNOWLEDGEMENTS

The research leading to the results described in this paper has received funding from the European Union 7th Framework Programme FP7-ENERGY-2009-1 under grant agreement number 241346.

The paper reflects only the authors' views and the European Union is not liable for any use that may be made of the information contained herein.