B Le-Roux

M Reimeringer

G Chantelauve

D Jamois

G Leprette

C Proust

L Mathieu

Industrial Risks and Land use Planning -Study of blast window resistance

Keywords: PPRT, window, blast loading test

Technological Risk Prevention Planning (PPRT) is a French tool for managing land-use planning near upper-tier SEVESO industrial facilities. Its purpose is to protect the population against industrial hazards. Risk limitation measures may include for example window reinforcement for new or existing blast risk exposed buildings. In order to better define technical guidelines for window reinforcement, INERIS developed an innovative experimental device. INERIS has studied window behaviour and has evaluated blast capacity of different kinds of glass panels

(monolithic or insulated glasses with or without anti-explosion film, laminated glasses...), mobile frames or window locking or anchorage systems. Results show that all window components have to be considered to improve safety.

Introduction

Most of experimental studies consider blast loading with high intensity and short duration simulating detonating devices. Moreover they have mainly focused on the behaviour of monolithic glass panels (annealed or tempered). Weisman et al. [START_REF] Weisman | Blast capacity evaluation of glass windows and aluminium window frames US Army Armament Research and Development Command[END_REF] or Giltaire and al. [START_REF] Giltaire | Essais de résistance de vitres aux explosions[END_REF], [START_REF] Giltaire | Comportement de vitres revêtues d'un film Scotchtint sous l'effet d'une onde de choc[END_REF], [START_REF] Mavrothalatassitis | Comportement des aux effets d'une onde de choc[END_REF] studied the response of monolithic or laminated glass subjected to blast shock wave within a short time (inferior to 10 ms). Coevert and al. [START_REF] Coevert | Blast tests on window panes treated with several types of experimental Shatter Resistant Film[END_REF] also studied security window film and insulating glass panel to high intensity blast load with a positive duration of approx. 1 ms to 10 ms.

In this paper, INERIS presented tests results of window subjected to a blast loading from accidental explosion characterized by a low intensity (20-50 mbar) and long duration (100 ms). These tests allowed to:

•

Evaluate blast load capacity of different insulating glass panels • Study the response of the whole window taking into account glass, frame, locking and anchorage system;

•

Identify what the weak points are;

Experimental device

Results analysis

The evaluation of glazing performance is done in accordance with the evaluate hazard rating criteria in Table I [START_REF]Glass in building -Explosion -resistant security glazing -Test and classification by shock-tube loading[END_REF] Table I Glazing shall be considered as "blast-resistant" only if it achieves a "minimal hazard" rating C or safer.

A complete window is considered as "blast-resistant" only if:

• Glass achieves a "minimal hazard" rating C or safer.

•

Window is still fixed to frame after the test • Mobile frame are still closed after the test • None piece of window is projected outside

Experimental results and discussion

Blast capacity resistance of glass panel

An initial series of [20-50 mbar] overpressure peak tests was performed to study the resistance of different glass panels. All panels were a 1.08 m x 0.60 m. Tests results are given in Table II and Table III. nnealed insulating glass panels resist an incident s of ms equivalent triangular duration. There w ncreasing overpressure to 35-50 mbar causes the ost fragments projected up to 3 m (Figure 9).

Blast capacity resistance of windows

A second test series studied the response of a French style inwards opening double window to a blast wave: resistance of the anchorage system, the fixed and mobile frame or locking system. All windows have a 1.25 m x 1.40 m 2 size. Test results are given in Table IV and Table V. The inward opening window brackets (Figure 13) which frame. Test results show tha overpressure blast wave. F window and a projection of For the second test the anc structural frame. Window resist either but the weak frames were broken and pr screwed to the building. T 15). 9 ndow has a standard anchorage system constitut hich are clipped to the window and screwed to that this system is not sufficient to resist a 20 Figure 14 shows deformation of angle brackets, on of the window between 1 m and 2 m. Given the results of previous tests, a French opening window with a reinforced locking system illustrated in Figure 6 was tested. The window is composed of wood frame, laminated insulating glass and a reinforced anchorage sysem. was subjected to a blast wave with a 50 mbar incident overpressure and an equivalent positive phase duration of hundreds milliseconds. Tests results showed that the configuration resists. The Glass panel was intact, The frame remained closed and attached to the frame (Figure 16).

Figure 16: French style inwards opening window with reinforced locking system, laminated glass and reinforced anchorage system subjected to blast wave -Photography after blast test

Conclusion

These experimental tests have produced new data on the response of windows subjected to a low incident overpressure and a long duration (hundreds of milliseconds) blast load.

Blast loading tests showed that a wet-glazed security window film improves the insulating glass panel resistance. Laminated insulating glass is also a good technical solution. However it is not sufficient to have a resistant glass panel. Reinforcing locking system and anchorage system are also recommended. All window elements have to be considered to improve safety.

Thanks to an experimental study and a theoretical analysis INERIS produced a practical guide [START_REF] Le-Roux | Guide pratique Fenêtres dans la zone des effets de surpression d'intensité 20-50 mbar, diagnostic et mesures de renforcement[END_REF] giving elements for improving the integrity of window subjected to a blast load. For instance, this guide gives recommendation on type and size of glass panel that can be used for various blast wave intensities and various blast wave types (deflagration or shock wave). It also indicates requirements for the locking system, the design and the number and position of the anchorage points.

Figure 1 :Figure 2 :

 12 Figure 1: Experimental device outsideface

Figure 3 : 1 F 2 2 Fig 3 resistance of 1 . 2 )

 3122312 Figure 3: Blast pressure gauge

  Figure 5: Standard locking system window

  ing

:Figure 7 :

 7 Figure 7: Cross-section through witness aera

ure 8 : 4

 84 Insulating glass subjected to a blast wave 3 -25 mbar was no damage the same glass to Insulating glass with a se internaly. Glass panel with During the blast The "Inne window and the outer pane (Figure10and Figure11). On the other hand the use In that case, the glass doe glazing hazard (B rate).

Figure 10 : 6 Figure 9 :

 1069 Figure 10: Annealed Insulati

Figure

  Figure 13: Window

  Figure 15: Projection of mofile frame of French style inwards opening window composed with PVC frame (at left) or wood frame (at right) and a common locking system (catch and strike)

Table II :

 II Glass panels subjected

		N°
		Glass
	Dimensions of glass pa
		Incid
	Characteristic of blast wave	peakoverp (mb Positive
		duratio
	Hazard level

Table III :

 III Glass panels subjecte

		cted to blast load (part 1/2) -Results of blast tests
			1	2
					4/1
			4/16/4	Annealed In
			Annealed Insulating glass	+ Wet-gla
					windo
		panel	h=1.08 m l= 0.60 m	h=1.08 m l= 0.60 m	h=1 l= 0
		cident
		erpressure	20-25	35-50	45
		bar)
		ive phase tion (ms)	> 500 ms	> 500 ms	> 50
			A	F
		cted to blast load (part 2/2) -Results of blast tests
		N°	3
		Glass	4/16/4 Annealed Insulating glass + daylight security window film	44.2 Laminate glass (ann
	Dimensions of glass pa panel	h=1.08 m l= 0.60 m	h=1 l= 0
		Incid cident
	Characteristic of blast wave	peakoverp erpressure (mb bar) Positive duratio ive phase tion (ms)	50-55 mbar > 500 ms	65 > 5
	Hazard level	F
	Table II shows that anneal
	overpressure for hundreds
	to the glass. However incr
	break (Figure 8) with most
		Figur
			5

Table IV :

 IV Window subjected to blast load (part 1/2) -Results of blast tests

	Dimensions of window	h=1.25 m x l= 1.40 m
	Opening window	French style inwards opening double window
			Standard anchorage system :
	Anchorage system	angle brackets clipped to the window and fixed to the	Reinforced system with angle brackets anchorage
			structural framing
	Locking system	Standard locking system (catch + strike)	Standard locking system (catch + strike)
	Frame		PVC	PVC
	Glass panel		4/16/4	4/12/44.2
		Incident
	Blast wave	peakoverpressure (mbar) Positive phase duration (ms)	20-25 > 500 ms	55-60 > 500 ms
	Results		Deformation and rupture of anchorage system Rupture of frame Opening the window Projection of the window between 1 m and 2 m	Opening the window Rupture of mobile frame Projection of the window casement between 2 m and 3 m

Table V :

 V Window subjected to blast load (part 2/2) -Results of blast tests

	Dimensions of window	h=1.25 m x l= 1.40 m
	Opening window	French style inwards opening double window
	Anchorage system	Reinforced anchorage system with angle brackets	Reinforced system with angle brackets anchorage
	Locking system	Standard locking system (catch + strike)	Reinforced locking system Individual closure of the openings
	Frame		Wood	Wood
	Glass panel		4/16/4	44.2-8-44.2
		60-65	65-70	65-70
	Blast
	wave	> 500 ms	> 500 ms	> 500 ms
	Results		Opening the window Rupture of mobile frame Projection of one of the opening at approximately 2 m	Mobile closed Glass panels intact frames	remain

Acknowledgments

The work presented in this paper has been supported by the French National Ecology Ministry (MEDDE), which is gratefully acknowledged.


J M Lacome

C Lemofack

J Reveillon

F X Demoulin

NUMERICAL MODELLING OF FLASHING LIQUID JETS DUE TO LEAKAGE OF LIQUEFIED GAS STORAGE

Keywords: Thermodynamic Equilibrium, flash-boiling, two-phase flow, droplet, vaporization

One of the missions of the INERIS is to assess accidental risks induced by industrial activities. The use of liquid gases is very widespread in industry but a large potential hazard exists in case of accidental release. These substances can be flammable and explosive like the LPG or toxic like anhydrous ammonia. In case of an accidental leakage due to a vessel breach or a line rupture, a violently change phase could take place inside theses products. This phenomenon called flashing jet are still bad understood and many strong hypothesis are often used to predict by modelling its physical consequences in the near field like thermodynamic behaviour of the whole jet, pool formation. However all these physical consequences have to be well evaluated in order to better estimate the explosive cloud or toxic cloud formation. The objective of the INERIS-CORIA work is to develop a new numerical model with the aim of simulating two phase jet resulting from a leakage in a pipe containing a liquefied gas. The work focuses especially on the effect of the vaporization and boiling process in the jet. A thermodynamic equilibrium model for vaporization was thus developed. To test the model, an atmospheric two-phase jet of butane, emanating from a circular orifice is considered. The modelling results show that the calculated temperature behaviour in the spray jet by comparison with the observations is generally satisfactory. This result cannot be obtained with classical vaporization model.

Introduction 1.Context

Accidental releases to the environment of pressurized liquid gases under ambient conditions could generate a large flow rate when the breach appears on the liquid filled containment. The superheated released liquid can form a two phase explosive or toxic cloud mixture. Potential consequences of these accidental releases are injuries, fatalities, destruction of installations and possible evacuation from the surrounding area of the accident. Examples of fatal accidents ( [START_REF] Board | Hazardous Materials Accident Summary Report Failure of Tank Car TEAX 3417 and Subsequent Release of Liquefied Petroleum Gas[END_REF], [START_REF] Board | Hazardous Materials Accident Brief[END_REF], [START_REF] Board | Railroad Accident Report, Derailment of Union Pacific Railroad[END_REF]) involving superheated liquids showed huge consequences. Predicting behaviour of the liquid-gas cloud mixtures by modelling is of direct relevance for industrial risk assessment.

At the breach level, in the case of an accidental release of pressurized liquefied gas, the product released in the ambient air is suddenly placed under temperature and pressure conditions that are such that a part of the liquid vaporizes violently. This phenomenon is generally called a "flash". In order to calculate safety perimeters around industrial installations, one of the objectives of atmospheric dispersion research projects of INERIS is to improve models of flashing releases in realistic industrial environments.

Equivalent source term models exist for flashing release in current long range dispersion models intended to predict toxic effects or explosive cloud formation. Several factors can, however, invalidate simplified equivalent source models, especially in the very near field where many complex phenomena can occur and where obstacles can be found.

We can easily say that total mass released and the release duration are major parameters. Others parameters as the velocity, thermodynamic state, and amount and droplet sizes of imbedded aerosols of the material at the exit of the rupture are also required as inputs to the jet and dispersion modelling.

A number of experimental, theoretical and numerical studies have been carried out [START_REF] Fauske | Source term consideration in connection with chemical accidents and vapour cloud modelling[END_REF], [START_REF] Boccardi | Two-phase flow through pressure safety valves. Experimental investigation and model prediction[END_REF], [START_REF] Wheatley | Discharge of liquid ammonia to moist atmospheres -Survey of experimental data and model for estimating initial conditions for dispersion calculations[END_REF] to perform calculations of two phase mass flow rates. Several empirical and analytical models were developed based on hypothesis related to saturation conditions. In order to have a better use of the models, they have been compared with measured data obtained during experimental tests performed by INERIS since the 90's [START_REF] Jérôme | Evaluation of two-phase flow models for accidental release and comparison with experimental data[END_REF]. These comparisons allowed us to conclude that to reach a good agreement between measured data and calculating flow rate is difficult but conservative models are clearly pointed out. This issue is relevant in a risk assessment.

Recent critical review [START_REF] Britter | Toxic industrial chemical (TIC) source emissions modeling for pressurized liquefied gases[END_REF] of term source modelling for toxic release scenario for pressurized liquefied gases showed there is "still a significant uncertainty in the overall modeling process". Main difficulties come from modelling suspension, evaporation of aerosol and rainout. The aim of the present work is to improve the acknowledge and the numerical modelling of aerosol evaporation.

Jets experiments

In order to improve model of flashing releases in realistic industrial environments, INERIS carried out two-phase butane and propane jet releases [START_REF] Bonnet | Experimental Study of Accidental Industrial LPG Releases: Rain out investigation[END_REF] into the atmosphere and ammonia [10] jet releases in a congested environment. Indeed most leakage problems from breakages of vessel or line rupture result in the formation of jets due to difference between the internal and the external pressure. These experiments involved to measure the main characteristics of the jet and to increase the understanding of the behaviour of superheated (flashing) liquid jets.

Modelling objective

Basic principles of flashing jets are summarized in Figure 1. The jet consists in 3 areas:

• The expansion zone (flash boiling, atomization): Here, the fluid expands from the vessel hole pressure to the atmospheric pressure. At the end of this zone, we assume that the jet consist only in gas phase and liquid phase. Both of them are at the boiling temperature. During the expansion, jet atomization occurs and the liquid jet ends here in droplets. • The entrainment zone (secondary break-up, droplets evaporation): the turbulent jet drives the ambient air. The energy brings by the air, of which the temperature is greater than that of the jet, is used for droplet vaporization in the two-phase jet. • The final dispersion zone: The entrainment by the air atmosphere heats the jet and decreases its velocity up to the wind speed.

The objective of this work is to develop a numerical model with the aim of simulating thermo dynamical behaviour of a two phase flow resulting from a breach or leakage in a vessel containing liquefied gas. A focus is made on the entrainment zone where a strong decrease of temperature was measured.

This paper aims to present the first results and comparisons between modelling calculations with experimental data.

A description of approaches used to simulate thermodynamic process inside flashing jets is given in Section 2. Results and discussion are presented in Section 3.

Methods

The first step of the work includes models used to describe the jet up to the end of expansion zone. These will be used as boundary conditions for the simulation of the jet from de beginning of the entrainment zone. In the second step, precisions are given on the new phenomena brought by the flashing liquid jet and the models suggested for it.

Expansion zone

To represent the flash boiling phenomenon that occurs at the exit of the injector up to the end of expansion zone, the isenthalpic Homogeneous Equilibrium Model (HEM) has been used. The model is used to determine the vapour mass fraction at the end of the zone.

The model supposes that the gas and liquid have the same velocity and the same saturation temperature due to the thermodynamic equilibrium assumption between the two phases. The mixture of the phases is homogeneous. These assumptions imply that:

• The two-phase flow is considered as a homogeneous fluid with properties at the middle of those of liquid and gas phase. • The flow quickly tends to a thermodynamic equilibrium (no energy exchange between the liquid-gas system and the outside). So it is supposed to have this thermodynamic equilibrium during the expansion. • Isenthalpic expansion: the enthalpy of the flow is constant while its entropy decreases.

There is another HEM model with Isentropic expansion, the entropy of the flow is constant while its enthalpy decrease. This causes the fall in temperature of the system [START_REF] Van Den Akker H.E.A | Discharges of pressurized liquefied gases through apertures and pipes[END_REF] showed that both of the assumptions are unrealistic but fewer errors are involved with the isenthalpic expansion.

In the case of this work, the HEM model presents the advantage of vapour/liquid mass fraction estimation. This will help to determine the mass flow of each phase from the total mass flow given by the experimental data.

Isenthalpic HEM Model:

( ) ( ) ( ) ( ) Cte = T h X + T h X = T h g l exp exp exp 1 * * - (1) ( ) dT C = T dh pl ∫ ∫ (2)
From the equation (1) describing the enthalpy conservation, we introduce (2) After a derivation by X:

( ) ( ) ( ) dX dT C C C X + X + dT C + dT C = T h T h pl pl pg pg pl i g i l * *         * - ∫ ∫ - - 1 (3)
By supposing X << 1 and C pg < C pl , we see that the term ( )

dX dT C X pl / 1 * * - - in (3) is equal to latent heat vaporization: ( ) ( ) dX dT C X = T Lv pl b * * - -1 (4)
The heat exchange due to vaporization is supposed greater than the heat due to warming:

Lv dT C pl < < *
. By taking account of this assumption after the integration of (4), we finally have:

Lv T T C = X i pl exp exp - (5)
At the end of expansion zone, the jet is in atmospheric pressure. So the liquid, which is in thermodynamic equilibrium with its vapour, is at its boiling temperature ( )

eb T = T exp .
Another issue of this work is the determination of the velocity induced by the pressure drop inside the injector and by the flash effect. The maximum kinetic energy can be estimated by the pressure drop and the variation of thermodynamic energy in the flow. Here, the velocity is considered to be driven mainly by the pressure drop, thus it is determined by the Bernoulli law. Additional velocity could be considered due to the flow expansion (liquid to gas) however there is also a reduction of velocity due to pressure loss.

All the above work is useful to estimate values in the jet at the end of expansion zone.

The simulations start at this point, since the CFD model used is the two-phase Euler-Lagrange description model. The experimental case presented in this paper comes from the data base of the FLIE (Flashing Liquids in Industial Environment) project Figure 2. Computation area [START_REF] Bonnet | Essais de rejets diphasiques à grande échelle[END_REF]. The numerical software used for this study is numerical package FIRE V8.41 from AVL.

The chosen case is whose of a vessel of liquid butane at 7.78 MPa. The liquid is released through a pipe and reaches the ambient air in a two-phase jet. This jet emanates from a circular orifice of 10 mm in diameter. The mass flow rate conservation between the orifice and the end of expansion zone give the diameter of this end zone : D = 60.45 mm. The numerical domain in which the butane is spreading is presented in Figure 2. It consists in a cylinder of 36 m in length and 6 m in diameter. To simulate a rejection of liquefied gas with wind effect, a co-flow of air surrounds the spray injection. The mass fraction of each phase is given by the HEM model. The boundary conditions are recapitulated in Table 1. Finally, the turbulence is modelled with the K-Ű model and the applied gravity is 9.8 m/s

Abramzon vaporization model

Vaporization is based on Abramzon vaporization model [START_REF] Abramzon | Droplet vaporization model for spray combustion calculations[END_REF]. By considering ρ and D the density, and the binary diffusion coefficient in the gas around the droplet, k and a, the droplet diameter and thermal conductivity, Sh and Nu the Sherwood and Nusselt number, the mass transfer rate m m is given by ( 6) and ( 7):

( ) [START_REF] Wheatley | Discharge of liquid ammonia to moist atmospheres -Survey of experimental data and model for estimating initial conditions for dispersion calculations[END_REF] Usually, the mass Spalding number M B (equation ( 8)) and the thermal Spalding number T B (equation ( 9)) are supposed to be the same ( )

1 ln + B πρDaSh = m M c (6) ( ) 1 ln + B Nu C ak π = m T c pg
B = B = B T M
. One part of the energy available in the gas phase goes for the evaporation of liquid phase. The other part goes for heating the droplet. This is described in equation (10) below.

s F, F, s F, M Y Y Y = B - -∞ 1 ( 8
) c pg s s T h C v ρ = B (9) ( ) 2 l s s s c πa Q + Lv v ρ = T T h - ∞ (10)
In the three equations above (( 8), ( 9), ( 10)), YF is the vapour mass fraction of droplet material in the gas phase, vs is the velocity in the gas phase and hc is the droplet convection coefficient. The subscripts s and ∞ mean the values at the droplet surface and far from droplet.

The introduction of the thermal Spalding number in (10) combined with equation [START_REF] Wheatley | Discharge of liquid ammonia to moist atmospheres -Survey of experimental data and model for estimating initial conditions for dispersion calculations[END_REF] give the heat quantity necessary for droplet heating. [START_REF] Abramzon | Droplet vaporization model for spray combustion calculations[END_REF] However thermodynamic conditions at the end of the expansion zone are very particular in the case of flashing jets. Thermodynamic equilibrium assumption imposes that the liquid and gas temperature are the same and equal to the boiling temperature. Additionally, near the inlet boundary conditions of the jet droplet, environment is composed of pure vapour ( )

( ) ( )         - - ∞ 1 ln + B Nu C Lv T T k a h πak = Q T c pg s c l
1 ≈ ≈ ∞ s F, F, Y Y
. Thus, this is a limit case between boiling and vaporization. As seen in equation ( 8), with these conditions, M B is uncertain. However, the security brought by putting the maximum value of the both vapour mass fraction at 0.99 is not sufficient. Thus, some simulation cases have shown that the classic Abramzon vaporization model is not suitable for post flashing two-phase jet. So there is a necessity to develop another model for vaporization. We propose a special procedure to solve this problem. It is considered a thermodynamic system (see Figure 3) at its initial state with certain liquid mass fraction ( ) • Mass conservation (12).

New model: Thermodynamic Equilibrium Model (TEM)

• Isolated system with constant pressure transformation: enthalpy conservation [START_REF] Bonnet | Essais de rejets diphasiques à grande échelle[END_REF].

1 = y = y k final k k initial k ∑ ∑ (12) ( ) ( ) Cte = y T H = y T H = H final k final k k k initial k initial k k k t * * ∑ ∑ (13) ( ) ( ) ( ) final final initial g initial l initial final T H T , T H = T f - (14)
The first step of this method is to find the enthalpy at the initial state with de quantities y and T in [START_REF] Bonnet | Essais de rejets diphasiques à grande échelle[END_REF]. Knowing that the enthalpy is constant, the second step consists in finding the new values of these quantities by applying a method for finding the zero of the function ( ) final T f (14). At the final state, thermodynamic equilibrium implies that the liquid and gas are at the same temperature. This method can be more interesting than the Abramzon model, because not only it resolve the particular problem of flash-boiling, but also it integrates the condensation phenomenon, which is supposed to occurs during the jet spreading.

The test of this model is done with the butane properties. The purpose of the test is to see if it is consistent by varying the initial quantities of y and T . In this paper we show the test of the model by varying the initial liquid mass fraction, with the liquid and gas temperature at 272.5°K and 350°K respectively. The total pressure is equal to 0.1MPa.

The Figure 4 shows the final temperature (left-aligned) and the final mass fractions of liquid and vapour (right-aligned) in the system. The model is consistent with the physics. If the liquid initial amount is small, because of the great temperature in the gas compared to the liquid temperature, it evaporates completely. The temperature evolution in this case of total evaporation is linear, accordingly to the ideal gas law. But, the system come to a stage where the energy in the gas phase is not sufficient to completely evaporates the liquid. Thus, as we can see in figure 5, which compare the vapour saturation pressure of butane at the final temperature and the butane vapour pressure in the system, there is an equilibrium between the two phases when liquid remain in the system. The next step is to introduce this model in the FIRE software. However, this model did not give information on how long it takes to evaporate. So, the model is completed to finally have the mass transfer rate between liquid and gas phase. This is done by adding a film of evaporation (Figure 6). It consists in taking an amount l β of liquid (15) and g β gas (16) which will be applied to the model. At the current state of the work, these amounts of liquid and gas are not yet determined precisely and have the same value β . In this paper, we will just compare the results given by taking account two different values of β .

Results and Discussion

It is complicated to have experimental results on the types of jets described above. But the INERIS were able to have more reliable measurements in temperature. These measurements show the same behaviour as seen on Figure 7, showing the evolution of the temperature in the jet axis with β=10% . The spray is surrounded by a hot gas environment. However, the temperature of the spray decreases up to a certain distance. Due to the initial droplet boiling temperature, close to the injection, the vaporization process dominates the flow. Since evaporation is an endothermic phenomenon, the spray jet cools down until there is no droplet enough. Thus the spray temperature rises only once the liquid vaporization does not have enough influence in the flow.

INERIS has done temperatures measurements with thermocouples at different points on three axis (up, horizontal and down) in the jet as seen on Figure 8. The comparison of these measurements with computation results (β=10%) in Figure 9 shows that the simulations are not far from reality. The differences can be explained by the fact that the choice of β and the calculations on the boundary conditions are not yet accurate.

Conclusion

The modelling results with the TEM model show that the calculated temperature behaviour in the spray jet by comparison with the observations is generally satisfactory. This result cannot be obtained with classical vaporization model.

As a perspective to this work, a model is under development to calculate the value of β by taking account the turbulence in gas phase, the time step, the droplet area.

Finally, the assumptions used to calculate the boundary conditions bring many lack of precisions. It would be desirable to do a CFD model of the material flow from inside the vessel to the end of expansion zone.

Figure 1 .

 1 Figure 1. Schematic representation of the areas of the flashing jet

Figure 3 :

 3 Figure 3 : Initial and final state in the thermodynamic system

  presence of other gas like air ( ) g y . The liquid temperature ( ) l T differs from vapour temperature ( ) g T . Like in HEM model, we consider thermodynamic equilibrium model at the final state. The physics assumptions here are: • Ideal gas.

Figure 4 .Figure 5 :Figure 6 :

 456 Figure 4. Test of the model of mass and enthalpy conservation by varying the initial mass fraction

Figure 7 :

 7 Figure 7 : Temperature evolution in the jet axis

Figure 8 :Figure 9 :

 89 Figure 8 : Temperature field in a vertical plane passing through the simulated jet axis with the experimental thermocouples positions : Th1-Th6 = up axis; Th7-Th12 = horizontal axis; Th13-Th18 = down axis

Table 1 .

 1 2 . Boundary layer conditions

	Variable	Inlet	Coflow	Outflow
	Material	Butane	air (gas)
		(liquid
		& gas)
	Velocity (m/s)	29.49	1	-
	Temperature (°K)	272.6	300.15	-
	Mass flow rate (kg/s)	1.33	-	-
	Vapour mass fraction	0.17	-	-
	Droplets size (µm)	100	-	-
	Turbulence intensity	10	1	-
	Length scale	0.005	0.1	-


Benoit Charrasse
email: benoit.charrasse@ineris.fr

Pierre Hennebert

Céline Tixier

Sediment-Water Partition Coefficients of Hydrophobic Persistent Organic Pollutants (POPs) for Different Sediments with Passive Sampler

Partition coefficients of polycyclic aromatic hydrocarbons (PAHs) and polychlorinated biphenyls (PCBs) were determined in sediments from five distinct contaminated sites in France (marine harbor, rivers canals and highway sedimentation tank). The assessment of the risk caused by such micropollutants requires, in most cases, the measurement of their availability. To access to this availability, low density polyethylene (LDPE) membrane samplers were exposed to these sediments under constant and lowlevel agitation for periods of 1, 2, 4, 7, 10, 14, 18, 22, 25, 29, 32, 39 and 46 days. Freely dissolved porewater contaminants concentrations were estimated from the concentration at equilibrium in the LDPE membrane. The depletion of contaminants in the sediments was monitored by the use of performance reference compounds (PRCs). Marked differences in freely dissolved PAH and PCB concentrations and resulting sediment-porewater partition coefficients between these five sediments were observed. Data set was tested onto different empirical and mechanistic models. As final findings, triple domain sorption (a total organic carbon, black carbon and oil phases model) could model PCB data successfully whereas the best fitting for PAH partitioning was obtained by Raoult's law model.

INTRODUCTION

Dredged contaminated sediment can often act as a sink for hydrophobic organic contamiants (HOCs) released into the surroundings environment. Strong sorption of nonpolar organic compounds to certain types of organic phase in sediments tends to reduce their potential of emission, their accessibility to microorganisms and their concentration in sediment porewaters [START_REF] Cornelissen | Extensive sorption of organic compounds to black carbon, coal, and kerogen in sediments and soils : Mechanisms and consequences for distribution, bioaccumulation, and biodegradation[END_REF]. This strong sorption, possibly to soot or char-like carbonaceous particles is responsible for the well documented failure of empirical relationships based on the fraction of organic carbon (f oc ) and organic carbon-water partition coefficients (K oc ) to correctly predict contaminant availability in sediments [START_REF] Bucheli | Quantification of the sootwater distribution coefficient of PAHs provides mechanistic basis for enhanced sorption observations[END_REF]. Furthermore, weathered mineral oil (MO) residues have been previously shown to exhibit sorption strengths of similar or even higher magnitudes as compared to black-carbon (BC), for polycyclic aromatic hydrocarbons (PAHs) and polychlorinated biphenyls (PCBs) [START_REF] Cornelissen | Extensive sorption of organic compounds to black carbon, coal, and kerogen in sediments and soils : Mechanisms and consequences for distribution, bioaccumulation, and biodegradation[END_REF]. Obviously, presence of soot and/or oil fractions in sediments may significantly limit HOC bioavailability and risks involved for the environment [START_REF] Cornelissen | Extensive sorption of organic compounds to black carbon, coal, and kerogen in sediments and soils : Mechanisms and consequences for distribution, bioaccumulation, and biodegradation[END_REF].

Our study aims to determine the equilibrium state of these chemicals for 5 dredged sediments representative of different types of contamination sources. The estimation of porewater concentrations (C pw ) of PAHs and PCBs were therefore determined by passive sampling experiments with low-density polyethylene (LDPE) strips exposed to the sediments in the laboratory. LDPE strips have been previously successfully used to assess freely dissolved concentrations of various HOCs in heavily contaminated harbor sediments [START_REF] Allan | Measuring nonpolar organic contaminant partitioning in three Norwegian sediments using polyethylene passive samplers[END_REF]. The depletion of the sediment by the passive samplers was monitored by the use of Performance Reference Compounds (PRC) (Brandt et al., 2012). Sediment-porewater partition coefficients, determined by combining these data with the total sediment concentration (C sed ), were then evaluated against various empirical and mechanistic partitioning models.

MATERIALS AND METHODS

Sediments Characteristics. The 5 dredged sediments were collected from contaminated sites across France and stored in 60 L barrels. All the sediment were selected to present contaminant concentrations higher than S1 threshold (for fluvial sediment) and N2 threshold (for marine sediment), according to French regulations on sediment dredging (Decree of August, 9, 2006). About 1200 kg (wet weight) of sediment were collected per site.

LDPE Membrane Preparation.

Lay-flat LDPE tubing was purchased from Brentwood Plastics Inc. (St. Louis, USA) and reproducibly cut along the two edges resulting in 2.5 cm wide LDPE membranes (average membrane thickness of 80 μm). 5 cm long membranes were pre-extracted twice with hexane, prior to the spiking with PRCs [START_REF] Booij | Spiking of performance reference compounds in low density polyethylene and silicone passive water samplers[END_REF]. The 10 PRCs were selected to cover a logK ow range from 4.53 (anthracene-d10) to 7.3 (CB 204). Strips were kept at -20°C until further treatment. LDPE Exposure in Sediment Slurries. Homogenized composites sub-samples were used for laboratory sediment exposition. Seventeen glass jars were filled with 350 g of wet sediment and one LDPE strip of 5 cm. These jars were covered with Teflon strip then capped with aluminium-lined lids and tumbled on a rotating system (5 rpm) over a period of 46 days at 20°C in the dark. LDPE strips were regularly sampled for kinetics study. Unexposed strips were also used as control samples over the exposure time period. C sed and C pw were both determined on the same samples.

Extraction and Analyses of LDPE Membranes and Sediment characterization. After exposure, strips were washed with milli-ro water and extracted twice by soaking overnight in hexane. Hexane extracts were combined and internal standards added. About 3 g of dry weight sediment sample were extracted using an accelerated solvent extractor (ASE 350,Dionex). Total organic carbon (TOC) concentrations were measured using a TOC-V analyzer (Shimadzu) coupled with a SSM-5000A. Black carbon was measured using the CTO375 chemothermal oxidation method [START_REF] Bucheli | Quantification of the sootwater distribution coefficient of PAHs provides mechanistic basis for enhanced sorption observations[END_REF]. Total hydrocarbons contents (THC) were determined using GC-FID (Autosystem XL, Perkin Elmer). PAHs and PCBs were analyzed by gas chromatography (Agilent 6850) coupled to tandem MS/MS mass spectrometer (Agilent 7000).

Overview of initial characterization of the 5 tested sediments is given in Table 1. Nonlinear and Linear Regression: Nonlinear and linear regressions were performed with Sigmaplot version 11 (Scientific Computing, USA).

RESULTS AND DISCUSSION

Contaminant Accumulation. A measurable accumulation for PAHs (14 among the 16 EPA) and PCBs (22 to 36 including the 7 indicator PCBs) into LDPE membranes was observed for all 5 sediments. Figure 1 presents the uptake curves obtained for 4 PAHs over the exposure time of 46 days. Accumulation process in the membrane obeys a first order kinetics for all the studied compounds and, as previously described (Huckins, 2006), the concentration in the membrane after a time t (C LDPE (t)) is expressed in Equation 1:

C LDPE (t) = C eq LDPE [1 -e (-k e . t ) ] (1)

where C eq LDPE (ng/g) is the concentration in the LDPE strip at equilibrium and k e (d -1 ) is the exchange rate constant. Data from recovery down to 85% and up to 115% were ruled out. PAH and PCB concentrations from control LDPE sampler were below limits of detection.

The uptake curves fit correctly with Equation 1 with significant coefficients of determination for all the contaminants. Triplicate analyses at day 14 and 46 show low relative standard deviation and good reproducibility (coefficient of variation between 1% and 9% for PAHs and 1% to 11% for PCBs). Most of the contaminants (92% and 96% of PAHs and PCBs respectively) had reached equilibrium state at 46 days. The exchange rates between LDPE strips and water were calculated from equation 1 for each contaminant. As shown in Figure 2, the estimated exchange rate constant k e is well correlated to the hydrophobicity of the compounds represented here by the logarithm of the octanol-water partition coefficient (logK ow ). The exchange rate constants for the added PRCs are in agreement with the values of the exchange rate constants observed for the native contaminants. Similar results were obtained by [START_REF] Booij | Dissolved PCBs, PAHs, and HCB in porewaters and Overlying Waters of Contaminated Harbor Sediments[END_REF] and [START_REF] Allan | Measuring nonpolar organic contaminant partitioning in three Norwegian sediments using polyethylene passive samplers[END_REF] with LDPE strips exposed to harbors contaminated sediment. The slight difference observed for sorption equilibrium between both contaminant families could be explained, for identical logK ow , by higher exchange rate constant values for PCBs (Figure 2). Moreover, as illustrated on Figure 3 for the PRCs, the exchange rates between LDPE strips and water differed for the five sediments tested despite the use of identical exposure conditions. This may be related to sediment properties such as particle size distribution, porosity, total organic carbon, black carbon and hydrocarbons content. Additionally, the level and origin of the contamination may influence the mass transfer into samplers.

Depletion Monitoring Using Performance Reference Compounds. Deducing porewater concentration from this kind of experiments requires to check whether a nondepletive situation was really obtained, i.e., whether the sorption capacity of the sampler was well below that of the sediment. The depletion was monitored by the use of PRCs as proposed by [START_REF] Brand | Advice on implementing bioavailability in the Dutch soil policy framework[END_REF]. A distribution factor (DF) corresponding to the ratio of the amount of PRCs remaining in the membrane and the amount of PRCs in the sediment at equilibrium was calculated for all the PRCs in each exposure conditions. When DF is 0.1 or less, the capacity of the coating is very small compared to that of the sediment. Consequently, the measured dissolved concentration in the porewater (C pw ) will not be affected significantly by depletion of analyte from the solid phase in the sediment. At the opposite, high DF values mean that there is little sorption capacity in the sediment and it is likely that the coating has taken a substantial fraction of the analyte from the solid phase of the sediment and therefore that the calculated C pw value does not represent the original value. We calculated the DF values for all contaminants from experimental PRCs DF values. We decided to conserve data when equilibrium state was achieved and when the DF values was 0.1 or less for PAHs and 0.3 or less for PCBs. An underestimation of the porewater concentration by a maximum factor of 1.3 could be thus expected for PCBs. According to our consideration of acceptance or exclusion data, data of PAHs with logK ow < 6.75 were conserved for the five sediments whereas data of PCBs with logK ow < 7.11 were conserved for all sediments except for Evry. 62% of DF values for Evry sediment exceeded 0.3. This could be due to a poor concentration of native PCBs in the sediment and potentially related to a sandy structure with a low rate of condensed organic carbon phase (e.g., Table 1). Partition coefficients K d and K TOC . We compared our observations of sedimentporewater partitioning to equilibrium expectations based on sediment organic carbon content (f oc ), assuming the well documented relationship K d =f oc *K TOC . We used logK TOC,model for PAHs and PCBs on the basis of the correlation with logK ow established previously [START_REF] Schwarzenbach | Environmental Organic Chemistry[END_REF]. The logK d,exp values observed were always higher than the theoretical values. The divergence of logK d,exp values from such predictions was highest for PAHs than PCBs and exceeded on average f oc K TOC products by about 1.5-1.9 log units. The K d,exp values for PCBs were higher than predicted by about 1 log unit. Absorption within organic carbon (f oc K TOC ) is not sufficient to explain the sorption of native organic contaminants in these five sediments, as previously reported [START_REF] Cornelissen | Extensive sorption of organic compounds to black carbon, coal, and kerogen in sediments and soils : Mechanisms and consequences for distribution, bioaccumulation, and biodegradation[END_REF]. These results highlight the need to consider other sedimentary phases to account for the observed, stronger retention of the compounds investigated. Nowadays, two kind of partitioning model are used in the literature: (1) mechanistic approaches mainly by accounting for the contribution of all sorbents, and (2) empirical approaches mainly by accounting for inner properties of each contaminants and sorption properties of one specific selected phase (coal tar, activated carbon) [START_REF] Arp | Estimating the in situ sedimentporewater distribution of PAHs and chlorinated aromatic hydrocarbons in anthropogenic impacted sediments[END_REF]. For mechanistic approaches, the in situ distribution of POPs was evaluated using the double and triple sorption domain model (Equations 2 and 3).

Kd = C s /C pw =f oc . K TOC + f BC . K BC . C pw n-1 (2) Kd = f oc . K TOC + f oil . K oil .+ f BC . K BC . C pw n-1 (3)
in which K d (L/kg) is the sediment-porewater distribution coefficient, C s (μg/kg) and C pw (μg/L) were the POP concentrations measured in sediment and water respectively, f oc , f oil , and f BC were the measured organic matter, oil, and black carbon weight fractions (Table 1). K TOC (L/kg), K oil (L/kg), and K BC ([μg/kg]/[ μg/L]) were sorption constants for these respective phases. These latter partition coefficients were calculated from compoundclass specific regressions with logK ow . Freundlich exponent (n) for nonlinear sorption to BC was set to 0.7 [START_REF] Hawthorne | Measured partition coefficients for parent and akyl polycyclic aromatic hydrocarbons in 114 historically contaminated sediments : Part 2. testing the kocKbc Two carbon type-Model[END_REF].

Recent relationship published between (logK BC or logK oil ) and logK ow are summarized in Table 2 (Equations 4 to 9). Experimental data from the 5 studied sediments were tested with double and triple domain sorption model, for the validation of the best combination among equation 4 to 9. Additionally, logK TOC value were calculated according to [START_REF] Schwarzenbach | Environmental Organic Chemistry[END_REF].

TABLE 2. Relationships between K BC , K oil and K ow for PAHs and PCBs.

For empirical approach, we tested the "coal tar model", considered as the most appropriate polyparameter linear free-energy relationships (PP-LFER) based on the Abraham's solvatation parameters [START_REF] Arp | Estimating the in situ sedimentporewater distribution of PAHs and chlorinated aromatic hydrocarbons in anthropogenic impacted sediments[END_REF]. Additionally, we also compared our set of data to the Raoult's model, using coal tar as sorbent phase For PCBs, and for both models, we added a specific term based on abundance ratio of PCB congeners (WF), as a proxy for weathering, aiming for a better correlation [START_REF] Hawthorne | Improving predictability of sediment porewater partitioning models using trends observed with PCBcontaminated field sediments[END_REF]. WF is ranging from 0 to near 1, respectively for fresh to weathered contamination. From experimental concentrations of CB 28, 52 and 153, we obtained WF of 0.60, 0.75, 0.83, 0.86 and 0.96, respectively for Evry, NB, Lens, Dunkerque and Marseille sediments.

Models were considered successful for the prediction of average logK d or logK TOC when values were at 90% fitted within a difference of ±1 log (factor ±10) and at 99% fitted within a difference of ± 1.5 log (factor of ±30) [START_REF] Arp | Estimating the in situ sedimentporewater distribution of PAHs and chlorinated aromatic hydrocarbons in anthropogenic impacted sediments[END_REF][START_REF] Hawthorne | Improving predictability of sediment porewater partitioning models using trends observed with PCBcontaminated field sediments[END_REF].

For PAHs, among tested models, the best results were obtained for empirical models and particularly for Raoult's law (Figure 4). This model predicts logK TOC values in the differences intervals of ±1 and ±1.5 with frequencies of 91% and 100% of the values, respectively. The deviation on average is -0.39 (±0.52). Finally, Raoult's law approach underestimates softly logK TOC predicted. On the other hand, mechanistic model exhibited better results for PCBs (Eq. 3 with estimation of K BC from Eq. 8 and K oil form Eq. 9). Observed values were all included in the factors 10 and 30 confidence interval (Figure 4). This approach presented high level correlation with experimental values, the deviation on average was 0,05 ± 0.32 (factor of 1.1).

CONCLUSION

Passive samplers with LDPE strips are suitable to measure porewater concentrations and thus to assess the sorption capacity (sediment/ porewater partition coefficient) of various contaminated sediments. Operational correlations between model predictions and experimental data were obtained for five distinct sediments. Predictive relationships for partitioning PAHs by mechanistic approach were proved to be insufficient. Best fit were obtained by Raoult's model using inner properties of each contaminant (water solubility) and coal tar sorbent phase concentration. For PCBs, triple domain sorption (a total organic carbon, black carbon and oil phases model) could model data successfully The differences observed in PCB partitioning in the 5 sediments could be explained by the type and source of organic matter and source of the contaminants. Nevertheless, considering these two models, it should be kept in mind that 86% and 53% of the PCBs and PAHs respectively, fitted within a factor of 3. In term of contaminated sediment management, despite a good prediction for PCBs, the availability of contaminants should be assessed in most cases.

FIGURE

  FIGURE 1. Accumulation of 4 PAHs (phenanthrene (), pyrene (), chrysene (), benz [a]anthracene () in LDPE membrane exposed to Lens sediment.

FIGURE

  FIGURE 2. Exchange rate coefficients (k e ) as a function of logK ow for passive sampling experiments with LDPE strips exposed to Lens sediment with PAHs (), PAH-PRCs (), PCBs () and PCB-PRCs ().

FIGURE 3 .

 3 FIGURE 3. First-order dissipation rate constant (log K e ) for PRCs spiked into LDPE membrane samplers exposed to Lens (), Evry (), Marseille (), NB () and Dunkerque () sediments.

FIGURE 4 .

 4 FIGURE 4. Representation of Raoult's law model for PAHs and triple sorption model for PCBs in studied sediments. The solid line is 1:1. The internal dashed lines represent the confidence interval suitable for a factor of 10 (1 log difference). The external dashed lines correspond to the confidence interval suitable for a factor of 30 (≈ 1.5 log difference).

TABLE 1 . Initial characterization of the 5 tested sediments including: black carbon (BC) and total organic carbon (TOC) content; contaminant concentration; total hydrocarbons content (THC) and texture.

 1 * naphthalene, acenaphthene, acenaphthylene, fluorene, phenanthrene, anthracene, fluoranthene, pyrene, chrysene, benz[a]anthracene, benz[a]pyrene, indeno[1,2,3-c,d]pyrene, dibenz[a,h]anthracene, benz[g,h,i]perylene *** CB 28, CB 52, CB 101, CB 118, CB 138, CB 153, CB 180

		Unit	Evry	Dunkerque	Lens	Marseille	Nimy Blaton
	Origin		HST*	Marine harbor	Canal	Marine harbor	Canal
	BC -black carbon	%	0.00	0.06	0.66	0.22	0.35
	TOC	%	3.20	2.85	13.9	4.93	10.6
	∑ 14 HAP**	µg/kg	2550	7710	45700	37100	3360
	∑ 7 PCB***	µg/kg	62.6	737	539	4050	312
	THC (C10-C40)	mg/kg	4590	3170	9580	5340	6310
	Texture		Sandy-silty	Silty-sandy	Clayey silt	Silty-sandy	Clayey silt
	*highway sedimentation tank

*

  Consequently, we only considered data of PCBs with logK ow values below 6.38 for this sediment. Our measurement of different PRC dissipation rates for different sediments exposed under identical batch experimental conditions indicate the importance of the PRC use to confirm the equilibrium state and avoid an underestimation of C pw .Freely Dissolved Porewater Concentration. C pw were then determined from LDPE concentrations at equilibrium fitted to Equation 1. For the present work, logK LDPE-water for both PAHs and PCBs were calculated from the molecular weight (MW) according to[START_REF] Smedes | Polymer-Water Partition Coefficients of Hydrophobic Compounds for Passive Sampling: Application of Cosolvent Models for Validation[END_REF] and considering, for PCBs, a specific term taking account of the chlorine position. The determined PAH dissolved concentration levels in porewater ranged from 0.146 ng/L for benzo[a]pyrene to 313 ng/L for pyrene. The dissolved concentration levels of PCBs were lower, ranging from 0.002 ng /L to 7.76 ng/L, respectively for CB 174 and CB 52. Higher PCB and PAH concentrations were observed, respectively for Marseille and Lens sediments.

ACKNOWLEDGMENTS

We acknowledge financial support from GEDSET INTERREG and ECOREV, and PhD granting from the Provence-Alpes-Côte d'Azur Region (France).


Jean-Marc Lacome

Benjamintruchot

HARMONIZATION OF PRACTICES FOR ATMOSPHERIC DISPERSION MODELLING WITHIN THE FRAMEWORK OF RISK ASSESSMENT

Keywords: safety distances regulatory studies, CFD, mass consistent, Gaussian model, integral model, blind tests, tridimensional (3D) models, harmonization

Since many years, using tridimensional (3D) models is increasing in order to compute hazardous consequences of dangerous phenomena such as fire, explosion or atmospheric gas dispersion and then estimate safety distances. In the specific context of the French regulation and mainly for the land use planning in the neighbourhood of industrial facilities as depicted in the Technological risk prevention plan (PPRT), prediction of safety distances is a key issue. First computations that were achieved has highlighted discrepancies between atmospheric 3D models (CFD, mass consistent, ...) results but also with conventional approach as Gaussian or integral models. In order to prevent such discrepancies and to enlighten scientific reasons of those differences, a French working group was created in November 2009 on the demand of French Ministry of Ecology and Sustainable Development. To reach these objectives, several simulation benchmark tests were carried out. In a first part we present the benchmark tests as well as their interpretations and we describe efforts of harmonization on practice and input data. In a second part we focus on the scientific analysis to perform more homogeneous input data between the several kinds of 3D atmospheric dispersion models. Some of the best practices and requirements are presented and discussed in the last part.

INTRODUCTION Context: Industrial risk management in France

On September 21, 2001, a huge explosion occurred in Toulouse (at the AZF factory) causing 31 deaths and thousands (~ 2500) of injured people. Two years after this major industrial accident, a new law was introduced on July 30, 2003 which described both prevention and repair of the damage caused by industrial and natural disasters. This law was the starting point of a general regulatory environment, then regulations have been made considerably tighter and the entire approach towards risk assessment has changed [START_REF] Lenoble | Introduction of frequency in France following the AZF accident[END_REF]. The Technological Risk Prevention plan (PPRT in French, standing for Plan de Prévention des Risques Technologiques) is the new legislation tool aiming to protect people by acting on the existing urbanization and also by controlling the future land-use planning in the vicinity of the existing Seveso establishments. The identification of potential scenarios is the first step of the numerous one required along the process of the PPRT elaboration. For a physically possible scenario, this analysis is then followed by the prediction of potential consequences (dangerous phenomena). The prediction of the impact area (thermal, overpressure and toxic effects), which is generated by a dangerous phenomenon, is another required stages. The prediction of effects distances is always performed using computational tools. Atmospheric dispersion models are then required to predict distances or impacted zones by toxic or flammable products. When the PPRT is finalized, it delineates a risk exposure perimeter, at the heart of which regulated zones are established. Among these zones, areas can be defined for expropriation or for protective measures on buildings which can be compulsory. Due to the new legislation the stakes are higher and this highlights the importance in computing precise distances to prevent from people exposure while assessing realistic safety cost. Predicting models used in this framework are characterized by a large variety of nature and complexity. For the same analysed hazardous phenomena considered, discrepancies appeared in terms of computed distances or impacted zones between atmospheric 3D models (CFD, mass consistent, ...; tridimensional model means other than integral or Gaussian models) results but also with conventional approaches as Gaussian or integral models. What are the parameters that can explain these discrepancies between those computations? Actually, reasons are numerous and of different types. The reasons generally put forward are given in the list below, however this list is not exhaustive: Differences between hypothesis and difficulties to be implemented as input data for dispersion models, Emission term source is estimated following distinct approaches, Intrinsic feature of dispersion atmospheric model : gaussian, integral, CFD, mass consistent, lagrangian, Regardless input data, implementation of atmospheric dispersion models of same feature may strongly differ, Results can be identical but interpretations in terms of effects distances are not equivalent.

Items 1, 3 and 4 are particularly quoted when 3D atmospheric dispersion models are used. Indeed, since many years, the use of 3D models is increasing in order to compute hazardous consequences of toxic/flammable gas dispersion. Particularly, one often argue that this kind of models allow to take into account complex environment (presence of obstacles) on the process of atmospheric dispersion in contrast with the traditional models. In order to prevent discrepancies between atmospheric dispersion modelling and to enlighten scientific reasons of these differences, a French working group was created in November 2009 on the demand of French Ministry of Ecology and Sustainable Development.

Objectives of the French tridimensional atmospheric dispersion working group

The objectives of this group were to explain scientific basis of 3D atmospheric dispersion modelling in order to justify links between results and input data or sensitive parameters, to create a guideline of best practices for 3D atmospheric dispersion modelling with the objective of forecasting hazardous consequences within the framework of risk assessment and to provide a reading tool for the administration in charge of the evaluation of the above mentioned studies.

METHODOLOGY AND PURPOSE OF TEST CASES

A 3D computation requires to carry out five main stages that are illustrated in Figure 1. In order to assess the influence of the whole set of calculation stages on final results, three tests cases, including two fully blind tests, were achieved by working group partners with their own model. Significant stages were identified: implementation of source term emission, choice of calculation domain, mesh generation, boundary conditions, choice of physical sub model. The term source release modelling is one of the most important steps in such an approach but not only for 3D modelling. However source term modelling was not considered within the scope of the study which focused on atmospheric dispersion modelling. The test cases computed by the different participants were used to elaborate good practices. About 12 modellers have participated to this exercise by using mainly k- based model. On top of that, other RANS and LES approaches and mass consistent models were also used. For each test, efforts of harmonization on practice and input data were performed. Those computations were used to highlight the main reasons of discrepancies and to elaborate new requirements

Sum up of test cases

As a first blind test, a fictive gas release in a free field was chosen. Three gases were studied: heavy, neutral and light gases. The set up of the release is expected to characterize a line rupture (2 inches diameter) on a vessel under pressure (8 bars). A ten minutes release was computed in the wind direction. Following French regulatory requirement (Circular of May 10, 2010) two meteorological conditions, which are characterized by Pasquill classification and wind module (F3 and D5), were modelled. As in regulatory studies, the definition of the first test case was limited to source term definition, atmospheric stability and wind definition. Modellers were free to set up others calculation parameters needed for the simulation. This first test case results clearly showed great discrepancies in the results, focused on toxic effects distances, because the way of using the 3D models was specific to each user. Some results are presented in Figure 2 and illustrate that for a low evolution close the toxic threshold a weak difference between calculations can lead to a substantial difference in terms of distance. This aspect is the basis of the hereafter discussion. It confirms a need to standardize the methodology, four major factors were identified: interpretation of input wind profile for 3D, turbulence models, mesh and cells size, source term implementation. The following blind test was performed in order to investigate the influence of these factors to prediction calculations. This second fictive test case consisted in modelling a massive release of propane (emission rate of 45 kg/s) in an obstructed field. This massive release is supposed to produce an unconfined vapour cloud explosion. A series of obstacles were introduced inside the domain. First, two factors were fixed: inflow wind, temperature and turbulence profiles, then, the source term was simplified. On top of global results comparisons, specific comparisons were made to evaluate the impact of some identified specific factor. Then, two main turbulence modelling approaches were evaluated: RANS and LES. In using similar code, differences on results were observed due to factors: buoyancy effects taken into account or not in the turbulence modelling, obstacles roughness modelling, surface or volume source term implementation, and the mesh. Although discrepancies were observed in terms of gas volume and explosive mass it is important to keep in mind that the final results associated to the explosive gas release, i.e. overpressure distances, are less sensitive to the dispersion modelling than prediction of toxic distances. Indeed, in some cases uncertainties on flammable cloud prediction calculations may sometimes be reduced by approaches used downstream to predict overpressure distances. This can be illustrated through the common approach, which is called Multi-Energy, to estimate overpressure distances that is based on a function of M ^1/3 , with M the total mass of explosive gas.

Specific works were performed during this second fictive test to investigate the consideration of turbulence production by buoyancy effects and distance upstream the first obstacles. Based on those two first blind cases, the working group observed that progress had been made by the different participants by using common practices. After having defined most of those good practices, a third case based on an experimental case (Kit Fox field campaign) has been carried out to evaluate it. A concise list of good practices is given in a later chapter.

STANDARDIZATION OF INPUT METEOROLOGICAL PROFILE: RELATION BETWEEN PASQUILL CLASSIFICATION PROFILE AND CFD APPROACH

As mentioned previously, for atmospheric dispersion modelling in case of ground level release, French regulation requires a minimal set of atmospheric conditions defined with the traditional Pasquill scheme: F3 and D5. A roughness length is also required to characterize the environment of the industrial plants. It is obvious that such conditions cannot be translated easily to a 3D model approach. Having in mind that inlet profile is one of the key parameter, it was then decided to set up some rules in order to achieve the homogenization of the inflow boundary conditions. The relation between Pasquill classes and tridimensional inflow boundary conditions is difficult to establish because the classes of Pasquill represent a broad variety of possible states of the atmosphere, as illustrated by the well known diagram of [START_REF] Golder | Relations among stability parameters in the surface layer[END_REF] in Figure 2. For an atmospheric condition defined only with three parameters (Pasquill class, velocity module at reference height, and ground roughness) several profiles are then possible. It is clear that, whether integral and Gaussian models are commonly using highly simplified velocity profiles, 3D models need more sophisticated profiles for different physical quantities as wind, temperature and turbulence profiles. It is also obvious that, in most of the industrial cases, not only the surface boundary layer must be modelled but the entire atmospheric boundary layer one. For regulatory reasons and consistence with existing computation, it is important to limit the set of input parameters meteorological data as much as possible. Great work was performed by the working group to establish a consensus between current theories [START_REF] Dyer | A review of flux-profile relationships[END_REF][START_REF] Panofsky | Atmospheric Turbulence[END_REF][START_REF] Zilitinkevich | Resistance and heat-transfer laws for stable and neutral planetary boundary layers: old theory advanced and re-evaluated[END_REF][START_REF] Gryning | On the extension of the wind profile over homogeneous terrain beyond the surface boundary layer[END_REF]. In previous guidelines intended for environment around buildings [START_REF] Tominaga | AIJ guidelines for practical applications of CFD to pedestrian wind environment around buildings[END_REF], formulations were suggested to express profiles in the surface boundary layer where the assumption of constant shear stresses is valid. From similarity theory, under horizontally homogeneous and steady state conditions the wind profile follows the logarithmic law (Panofsky and Dutton, 1984, p. 134):

                 L z z z u z u M 0 0 * ln ) (  (1)
Here κ (0.4) is von Karman constant, z is height above the ground, u *0 is friction velocity friction, u is mean wind speed, M are the universal function of the stability parameter z/L, where L is defined as the Monin Obhukov length. While it gives a valid description of the surface boundary layer, this formulation is rigorously not valid above the surface boundary layer. In the scientific literature there are fewer proposals to extended profiles above this latter. As indicated previously, the working group has to pay a specific attention to these extensions in order to simulate a more realistic boundary layer flow. The approach proposed by the working group is based on the Gryning theory [START_REF] Gryning | On the extension of the wind profile over homogeneous terrain beyond the surface boundary layer[END_REF] and could be described as follow:

1. First, following former approaches used for integral and Gaussian models, only 4 parameters are required: Pasquill Class, Wind module (u ref ) at height reference (z ref ), length roughness (z 0 ), ground temperature (T 0 ); 2. Next, the Golder approach allows to determine a value of Monin-Obukhov length. To establish a unique value of Monin-Obukhov length, it was decided to rely on formulation that express midpoint of Pasquill class:

) ( log 1 1 0 10 s S z z L L  (2)
where for range 0.001  z 0 .0.5, parameters Ls(m) et zs(m) are constants. 3. Finally, the Gryning approach is used to express wind profile within the surface boundary layer profile and above surface layer. This leads to an iterative calculation to estimate the friction velocity near the ground, u *0 that matches with required wind module at height reference; kinetic turbulent energy and dissipation profiles are based on the equilibrium hypothesis. In the surface layer, the friction velocity u *0 , is considered constant as previous authors theory [START_REF] Businger | Flux profile relationships in the atmospheric surface layer[END_REF][START_REF] Dyer | A review of flux-profile relationships[END_REF]. Above the surface layer the local friction velocity has been considered to diminish with height, the following expression is used:

) 1 ( ) ( 0 * * z z u z u i   (3)
with z i the height of the atmospheric boundary layer. In order to keep a minimal of turbulent energy close to the top of boundary layer it was decided to consider a minimal value for u * (z) equal to 10% of ground local friction velocity. The wind profile is expressed following the most general formulation of Gryning's approach:

)) 2

( ) 2 1 ( ) (ln( ) ( 0 0 * L z z z L z z z L bz z z u z u MBL i MBL i       (4)
with L MBL that express length scale in the middle of the boundary layer and z i the top boundary layer. For stable conditions the most general formulation of L MBL , can be expressed following formulation:

)] 2 ( ) ln( )] ( )) ( ) [[(ln( 2 0 2 / 1 0 * 2 2 0 * 0 0 * L z b z z fL u A fL u B fz u L z i i MBL i      (5)
with f the Coriolis parameter (f=1.10 -4 s -1 ). Gryning pointed out that these coefficient values are rather well established for the neutral atmosphere but their stability dependence is still matter of discussion [START_REF] Zilitinkevich | Resistance and heat-transfer laws for stable and neutral planetary boundary layers: old theory advanced and re-evaluated[END_REF]. Coefficient values of A = 10 and B = -10 were taken from [START_REF] Tennekes | Similarity laws and scale relations in planetary boundary layers[END_REF] by and correspond to stable cases. The working group concentrated its effort on this kind of conditions because they are conservative for a majority of regulatory studies but the whole approach can be adjusted to others conditions (neutral, instable). Above the top layer the wind velocity is fixed to value obtained at z= z i . This latter is estimated by an analytical formulation [START_REF] Garratt | The atmospheric boundary layer[END_REF]:

f L u zi 0 * 4 . 0  (6)
To evaluate the turbulent kinetic profile, the similarity functions proposed by [START_REF] Dyer | A review of flux-profile relationships[END_REF] for stable conditions were considered, the profile is specified as follows:

L z L z C z u z k / 5 1 / 4 1 ) ( ) ( *     (7)
The dissipation rate profile is based on hypothesis that viscous dissipation balances shear production and buoyancy. Hypothesis of local friction velocity was considered consistently with wind profile formulation.

) 4 1 ( ) ( ) ( 3 * L z z u z     (8)
The turbulent viscosity, K m , is expressed by: )

( ) ( ) ( 2 z z k C z Km    (9)
with C µ the coefficient use to define the eddy viscosity in the k-ε model. Standard value for C µ = 0.09 and value recommended by Dyunkerke [START_REF] Duynkerke | Application of the E-ε Turbulence Closure Model to the Neutral and Stable Atmospheric Boundary Layer[END_REF] were tested by modellers. An application of this approach is being tested thanks the third test case performed by working group by comparisons with an experimental case. The Kit Fox field campaign, trial 5-8, that corresponds to stable conditions was chosen. As explained previously a set of minimum data was used: u ref = 4.1 m/s, z ref = 4,88 m, z o = 0.2 mm, T 0 = 301.8 K, L = 8 m.

DISCUSSION ON BEST PRACTICES AND REQUIREMENTS

A list of numerous best practices, requirements on 3D dispersion model capabilities and users' skill to demonstrate it has been listed by the French working group. Some practices come from its own work: construction of the inlet velocity profile, definition of a validation procedure, range values of Schmidt coefficient and justification,...; others were already identified in previous guidelines [START_REF] Franke | Best Practice Guideline for the CFD simulation of flows in the urban environment[END_REF] as the minimum distance between obstacles and domain boundaries, etc.

CONCLUSION

In this paper we presented the harmonization of practices for atmospheric dispersion modelling within the framework of risk assessment. The different classes of stability, required by the regulation, have to be constructed on the basis of equations established by the working group. If we focus on suitability of 3D models to simulate atmospheric turbulence in order to be in accordance with regulation, the demonstration of model capability to maintain wind and turbulence profiles through the free field calculation domain, up 2 km limit, is one of the main requirements. Knowing the complexity for CFD model to fulfil this requirement in a strictly way, it was decided that outlet profiles must show the same Pasquill class as input profiles. To conclude, it is worth noting that for some chemical products many uncertainties still exist in terms of acute toxic threshold values and scientific researches are going on this topic, which 3D models could allow foreseeing improvements by taking into account the intermittency of toxic cloud.

Figure 1 :

 1 Figure 1 : Main calculation stages to simulate atmospheric dispersion of an accidental release (toxic or flammable gas release)

Figure 3 :

 3 Figure 3: Relations between Pasquill stability categories and Monin-Obukhov length and surface roughness (from[START_REF] Golder | Relations among stability parameters in the surface layer[END_REF]

  Figure 2 : calculated concentration of CO at 50 m, 100 m, 500 m for two calculations of the first test case ; SEI 10 min standsfor Seuil effet irreversible is the threshold of irreversible effects for 10 min exposure; SEL 1% stands for Seuil effet létal is the threshold of first lethal effects above which 1% mortality can be observed in the exposed population

	Term source	1	2	Grid mesh	3	4	Post processing
	release modeling : jet under pressure, pool	Implementation of term source release		design, boundary layer conditions	Atmospheric sub physical model		: toxic dose, explosive cloud, effect distances
	evaporation, ...

ACKNOWLEDGEMENTS

Authors acknowledge the overall French working group and are particularly grateful to members who contributed to the scientific activities: B. Carissimo, J. Commanay, L. Soulhac, A. Albergel, C. Souprayen, A. Stanissopoulos, A. Carrau, M. Ichard, Marc Till, A. Champassith.


J M Lacome

F Tognet

C Olry

M Nibart

NEW LAGRANGIAN APPROACH FOR WET PLUME MODELLING

Keywords: Microphysics, lagrangian, dispersion, modelling, spray, microorganisms, legionella, wet plume

A new version of the lagrangian dispersion modelling system Micro-SWIFT-SPRAY (MSS developed by ARIA Technologies and ARIANET) has been developed in order to allow simulation of wet plume, microorganisms, insecticide or pesticide dispersion. For all these applications, it is necessary to take into account both gas and liquid phase evolution including mass and heat transfer. In the particular case of microorganisms distinct behaviors potentially occur when transported by the liquid and/or the gas phase. Following the MSS lagrangian approach, it is easy to define a mass quantity of vapor and a large number of droplets of the same diameter by a virtual particle. Hence a spectrum of different size of droplets can be model using several virtual particles. The development of the new module consists of modelling the physical phenomenon of evaporation and condensation inside a wet plume taking into account temperature and humidity of ambient atmosphere. Microphysics of droplet is solved using the classical laws of evaporation/condensation processes in the surrounding atmosphere. The calculated evaporation rate allows to estimate the diameter evolution and the temperature inside the droplet. The difficulty consists in characterizing the surrounding atmosphere. The interaction between droplets and ambient humidity has been performed on an eulerian frame. Two methods of two way coupling have been tested on academic cases (0D). The proposed paper will present the two approaches and the obtained results.

INTRODUCTION

Microphysical processes acting in a wet plume can be taken into account in 3D dispersion models like CFD. However, simulations can be complex to carry out and costly in computation time, particularly when a release have to be simulated over a few days to several weeks. Several Gaussian models offer simplified microphysics modules, however these models include simplified description of physics which is not relevant for complex flows which occur in urban areas. Even so applications may be multiple: transport of bioaerosols, plumes visibility and influence of wet plumes on plume dispersion. The development of a microphysics module in the Lagrangian model Micro Swift Spray allows to overcome the lack of tools available for these issues.

The basic of the new model is that virtual particles created by the lagrangian model include droplet spectrum of liquid water, vapor and microorganisms. Initially, the microorganisms are contained in water droplets. Evaporation rate is calculated for each droplet using the number of sherwood and correlation of Ranz and Marshall (1952 a,b). The evolution of droplet temperature is driven by the generic Lagrangian equation [START_REF] Miller | Evaluation of equilibrium and non-equilibrium evaporation models for many-droplets gas liquid flow simulation[END_REF]. Water vapor diffusion into the air is solved on a eulerian grid, this allows computation of diffusivity using the formulas of List (1963) and [START_REF] Bolz | CRC Handbook of tables for applied Engineering Science[END_REF]. Finally, the evolution of the density and temperature within the plume is determined by using the equations of [START_REF] Glendening | Improving plume rise prediction accuracy for stable atmospheres with complex vertical structure[END_REF]. When Water droplets are completely evaporated, microorganisms are released and start to disperse as free airborne. A simple biological model based on an exponential survival of microorganisms is applied at this specific moment when microorganisms are directly exposed to radiation, air temperature and humidity.

In the specific case of Legionella, accidental spread in the atmosphere due to contaminated cooling towers (CT), may disperse over a large distance [START_REF] Nguyen | A community-wide outbreak of legionnaires disease linked to industrial cooling towers--how far can contaminated aerosols spread[END_REF]. In this work a specific biological module has been developed and coupled with the MSS microphysic module. A first evaluation is done with the simulation of dispersion field campaign of biological aerosols performed by CSTB, DGA and INERIS.

MICROPHYSIC OF DROPLET

A multiplicity of thermo dynamical processes can occur within a plume released by a cooling tower. At the tower exit the plume is composed by a turbulent mixture of liquid droplets and vapour gas, which temperature is generally greater than ambient temperature. This bi-phase mixture is generally emitted with a vertical velocity that ranges from ~ ms -1 to a dozen ms -1 for large towers and as the plume rise it will be transported by wind. The liquid fraction value at the tower exit could be approximately considered around ~ 10% or less. This fraction is composed by a spectrum of droplets which evolves due to phenomena of nucleation, coagulation, evaporation, condensation, all of them being influenced by the turbulence of the flow. The droplets can coagulate between them mainly within area of the plume with high number densities. Such a phenomenon contributes to decrease the number of droplets and to reach a diameter of equilibrium. Heat transfers between droplets and vapour contribute to the evaporation of the droplets by heat mass transfer depending of the saturation level. Far away the droplets are assumed to reach an equilibrium temperature with the ambient air. The kinetics of evaporation and condensation strongly depends on the humidity saturation of the surrounding atmosphere. This latter evolves following interactions between the wet plume and the ambient flow which humidity level is a variable of primary importance. At the emission exit the droplets which are large enough can be affected by gravity, particularly those located in the surrounding of the release. In case of presence of drops with large radius (> 40 µm) [START_REF] Bouzereau | Application of a semi-spectral cloud water parameterization to cooling tower plumes simulations[END_REF] precipitation phenomena could occur but it's unlikely the case in the investigated applications (CT). In the far field dispersion, the median diameter is smaller but sedimentation can also affect the vertical velocity of the droplet.

Due to the great number of phenomena described above it was decided to launch the development of the new model by focusing on modelling evaporation and condensation. These phenomena play major roles in the temperature and diameter equilibrium of the droplets and its lifetime in the atmospheric far field dispersion.

NEW LAGRANGIAN APPROACH FOR WET PLUME

In the LPD model Micro-SPRAY, the dispersion of a wet plume is simulated following the motion of a large number of fictitious particles, each one representing a part of the mass of vapour and liquid water emitted from sources of general shapes. Starting from a given spectrum of droplets, each virtual particle carries a large number of droplets sharing the same geometric and thermodynamic properties. Particle's transport is obtained applying an equation of motion where the particle velocity is split into a mean component driven by the local wind reconstructed by Micro-SWIFT and a stochastic component obtained by solving a 3-D form of the Langevin equation of the random velocity [START_REF] Thomson | Criteria for the selection of stochastic models of particle trajectories in turbulent flows[END_REF].

As far as the wet plume has a vertical or oblique initial momentum and an initial buoyancy, five governing conservation equations of mass, energy, vertical momentum and the horizontal momentum are integrated for each particle at each time step, based on [START_REF] Glendening | Improving plume rise prediction accuracy for stable atmospheres with complex vertical structure[END_REF], [START_REF] Hurley | The Air Pollution Model (TAPM) Version 3. Part1: Technical Description[END_REF] and [START_REF] Hurley | Plume Rise and Enhanced Dispersion in LADM[END_REF]. These equations are described as follow:

2 2 2 p p p s w v u u    (1)
with u p , v p, w p velocities of particles at emission point.

 

a p a g B      (2) e bu E 2  (3)                         2 1 2 2 2 1 cos cos 1 cos cos a p a p a a p a p a s e U U u u           (4)
with , p refer to air and plume, respectively, b is the plume radius (m), B is the buoyancy (ms -2 ), E represents the entrainment rate (m 2 s -1 ), u e is the entrainment velocity (ms -1 ) and Ψ and Φ are the vertical and horizontal angles. The following conservation equations are solved:

mass s s a p u E b u dt d        2   (5) energy   2 2 2 b w u N B b u dt d p s a p s     (6) vertical momentum s p s a p u b B b w u dt d 2 2          (7)
x horizontal momentum

a s p s a p u u E u b u dt d          2   (8) y horizontal momentum a s p s a p v u E v b u dt d          2   (9)
where u a and v a are the horizontal components of the wind velocity, N 2 is the Brunt-Vaisala frequency. On one hand, by integrating this set of conservation equations, the temperature of the plume at the location of each virtual particle is known and used to drive the evaporation process of droplets. On the other hand, the concentration of water vapour is computed on an eulerian grid at each time step, using the background 3D field of specific humidity provided by Micro-SWIFT [START_REF] Moussafir | A new operational approach to deal with dispersion around obstacles: the MSS (Micro-Swift-Spray) software suite[END_REF] and the vapour carried by fictitious particles. In each grid cell, concentration of vapor is compared to the saturation concentration and, as a result, evaporation or condensation can occur.

In grid cell where condensation or evaporation takes place, transfer of mass between vapor phase and liquid phase occurs internally to each virtual particle present in this grid cell. Where condensation happens, this transfer is done until the concentration of vapor is locally equal to the saturation concentration. Where evaporation happens, two different approaches can be followed: 1) Microphysic of droplets disabled: Transfer from liquid phase to vapor phase is done until the concentration of vapor is equal to the saturation concentration. 2) Microphysic of droplets enabled: evaporation rate at each time step is calculated for each type of droplets using the number of sherwood and correlation of Ranz andMarshall (1952 and1953) together with the molecular diffusivity. This diffusivity can be computed by using the formulas of List (1963) or [START_REF] Bolz | CRC Handbook of tables for applied Engineering Science[END_REF]. Given this evaporation rate, mass and diameters of droplets are updated. Finally, the evolution of droplets temperature is also updated by using the following equation [START_REF] Miller | Evaluation of equilibrium and non-equilibrium evaporation models for many-droplets gas liquid flow simulation[END_REF]:

m dt m d C L T T Nu dt T d droplet droplet L v droplet air d droplet ) / ( ) )( ( Pr 3 1      (10)
where dm droplet /dt is negative for evaporation, T air and T droplet are respectively the local carrier plume temperature and the droplet temperature, L v is the latent heat of evaporation, the ratio of the gas heat capacity to that of the liquid phase is 1  =C p,air /C L , Pr is the gas phase Prandtl number;

VALIDATION TEST (0D )

In order to validate the evaporation process when microphysic of droplets is enabled in MSS, the evaporation rates of several sizes of droplets subjected to different humidity, temperature and wind velocity conditions were computed. Results have been compared to measurements of evaporation rates and model predictions both described in [START_REF] Kincaid | A Water Droplet Evaporation and Temperature Model[END_REF]. His model is based on modified Ranz and Marshall equations proposed by [START_REF] Goering | Mathematical modelling of spray droplet deceleration and evaporation[END_REF]. Figure 1 shows that the results obtained when using List diffusivity formula are in good agreement with both Kincaid's model and measurements. Starting from initial droplet temperatures situated above and below the wet bulb temperature, the evolution in time of the temperature of the droplet was tracked until |dT droplet /dt| < ε (where ε is a small threshold value). Drops of temperature predicted by MSS and by the model used in [START_REF] Kincaid | A Water Droplet Evaporation and Temperature Model[END_REF] are shown on Figure 2. Results from MSS calculation show that the temperature evolution is very similar to the Kincaid model prediction. However we notice that in MSS the droplet temperature drops below the wet bulb temperature. Shirai et al. (Shirai et al., 1971a,b) emphasized that measured droplet temperatures slightly lower than wet bulb temperatures were experimentally observed. An explanation provided in [START_REF] Kincaid | A Water Droplet Evaporation and Temperature Model[END_REF] could be that a water droplet may evaporate more efficiently than a wick-covered thermometer.

APPLICATIONS FOR A BIOLOGICAL MODEL

The model described above has been extended to model microorganism's transport. The underlying principle is that the virtual particles created by the lagrangian model include a droplet spectrum of liquid water, vapour and microorganisms. Initially, the microorganisms are assumed to be contained in water droplets. When water droplets are totally evaporated, microorganisms are released and transported as free airborne tracers. With respect to microorganisms, a consideration of sustainability must also be included in the code. An experimental laboratory approach specific to Legionella pneumophila has been carried out by CSTB and provided an useful knowledge (Ha T.L, 2005) on the bacteria sustainability using radiation, air temperature and humidity parameters. Following this experimental approach a model of sustainability has been set and coupled with the transport model. This biological model based on an exponential survival of microorganisms is applied during the free airborne transport only.

The main characteristics of this biological model are as follow:  model is disabled for bacteria located in a virtual particle including water liquid,  model is enabled for bacteria located in a virtual particle as soon as this particle is totally dry. At any time step t  , if one or more of these conditions of sustainability is not satisfy, a number of spores lost by the virtual particle concerned is estimated. Naming i  the decay parameter of the law of sustainability, the number of spores remaining after application of the survival model is as follows: A study supported by the AFSSET (Agence Française de Sécurité Sanitaire de l'Environnment et du Travail) consisted of estimating the ability of the model Micro Swift Spray to predict the potential spread of microorganisms by an anthropogenic source resembling a cooling tower, by comparing data from the simulation with experimental measurements. Results were presented in the HARMO 14 session [START_REF] Tognet | Numerical modelling of microorganisms dispersion in urban area: Application to legionella[END_REF] and parts of the conclusion are presented here. The simulation of the health risk caused by L. pneumophila aerosolization from a CT is possible but would require a finer characterization of the environment, the air flow conditions and source term, as well as additional knowledge on the ecology of Legionella before and after aerosolization (isolated bacteria in biofilm or with other interaction), the latter may be critical to the sustainability model. The choice of the model type (i.e. Gaussian, Lagrangian or Eulerian, with specific sustainability model or not) will depend on the expected accuracy and duration of the simulations requested by the authorities for the decision making process. It was experimentally demonstrated that, the dispersion of relative efficiencies of biocollectors is three times higher for these tests compared to controlled environment. While this does not preclude the ability of biocollectors to detect low concentrations of aerosolized microorganisms, the development of metrology is an important goal for improvement.

DISCUSSION AND CONCLUSION

The purpose of the present work is to develop a dispersion model for better understand outbreak consequences from a suspicious source of wet plume release over a few days to several weeks. If we consider that the modelling system Micro-SWIFT-SPRAY ensures a reliable description of the meteorological process in a reasonable computational effort, a model of evaporation/condensation for droplet has been implemented within the lagrangian dispersion model SPRAY. The droplet temperature is driven by the generic equation described in [START_REF] Miller | Evaluation of equilibrium and non-equilibrium evaporation models for many-droplets gas liquid flow simulation[END_REF] In order to validate (0D), evaporation rates of droplets, subjected to different conditions, were computed and comparisons with both Kincaid's model and measurements [START_REF] Kincaid | A Water Droplet Evaporation and Temperature Model[END_REF] give good agreements. A procedure of test case in 3D is presently carried out. By comparisons with both other lagrangian approach implemented in SPRAY [START_REF] Mortarini | Development of the microscale lagrangian particle dispersion model MicroSpray for the simulation of two-phase releases[END_REF] and CFD Code_Saturne [START_REF] Bouzereau | Application of a semi-spectral cloud water parameterization to cooling tower plumes simulations[END_REF]. These comparisons should provide hindsight on both lagrangian approach and eulerian approach.

d

  the particle time constant for Stokes flow and Nu the Nusselt number.

Figure 1 :

 1 Figure 1: Evaporation rate of droplets in different conditions

Figure 2 :

 2 Figure 2: Evolution of temperature of droplets of different sizes, with initial temperature below or above the wet bulb temperature; left figure: prediction by MSS, right figure: Kincaid et al. (1989)



  When sustainability model is enabled, sustainability criteria are as follow: o Minimal and maximal temperature of sustainability (T min and T max ) o Net minimal and maximal radiation (Ray min and Ray max ) o Minimal and maximal relative humidity (HR min and HR max )For each -time step, the local temperature T, the local radiation Ray and relative humidity RH (taking into account the relative humidity provided by Meteorological input as well as the contribution of the water vapor contained in the wet plume) are checked with the following conditions :