Frédéric Antoine

Jérôme Hébrard
email: jerome.hebrard@ineris.fr

Jean-Marc Lacome

ASSESSMENT OF THE MODELS FOR THE ESTIMATION OF THE CO2 RELEASES TOXIC EFFECTS

Currently, INERIS is involved in European and French projects regarding the CCS chain (carbon capture and storage). Nowadays, some people consider this chain as a future device to fight against the global warming due to high concentration of CO2 in the atmosphere. However, in case of massive accidental leak of CO2, this substance could be the origin of toxic effects for human. Now, carbon dioxide concentrations considered as potentially toxic get close to 100 000 ppm. To estimate precisely the distances reached by this hazardous level of concentrations, an efficient understanding of the CO2 release phenomena, from the assessment of the mass flow rate to the atmospheric dispersion, is necessary. Whereas the carbon dioxide is often stored and handled under 2-phase or supercritical conditions associated to storage pressure, CO2 ice formation is possible in case of accidental leak in the atmosphere due to specific properties regarding its triple point. Then, this CO2 flakes creation may be followed by the liquid/solid CO2 pool formation on the ground. Due to the important pressure drop of the fluid during the leak, a significant expansion phase, a high rate of air entrainment and a huge temperature drop followed by a dense and cold cloud formation should be considered. Concerning these specific points, a few of atmospheric dispersion softwares take into account the carbon dioxide specific conditions of release. And moreover, some calculations carried out with other software seem to provide different results in terms of hazardous distances for instantaneous or continuous releases. A focus on the models used by the softwares appears interesting to let a better understanding of this difference in the risk assessment context.

INTRODUCTION

Among several topics concerning the industrial activities, INERIS is involved in quantification of the hazard in case of accidental toxic releases. Indeed, hazardous substances, like chlorine, ammonia but also acids, LPG, hydrocarbons and so on, may be used and handled on industrial facilities. In case of severe accident, it appears possible that a big amount of substance may be released in the plant but outside as well. In such a situation, it is interesting to assess as precisely as possible concentrations of the substance in the air at different distances to predict potential toxic effects on people located around the site. Since few years, another substance is focused on due to global warming considerations: carbon dioxide (CO2). Some people consider it is possible to reduce carbon dioxide concentration in the air by trapping it in huge quantities with adapted facilities. In this case, storage, handling and transport of this substance should be considered in terms of toxic hazards in case of massive release. Indeed, high level of CO2 concentrations could be the origin of toxic effects for human. Currently, carbon dioxide concentrations considered as potentially toxic get close to 100 000 ppm [START_REF] Harper | The inclusion of CO2 as a hazardous substance in the Seveso Directive[END_REF]. Weaker levels could be interesting as well. Currently, INERIS is involved in European and French projects regarding the CCS chain (Carbon Capture and Storage), and more particularly pipeline transport. One of the main goals of these projects is to compare theory and results from experiments in order to estimate precisely the distances reached by this hazardous level of concentrations. An efficient understanding of the CO2 release phenomena, from the assessment of the mass flow rate to the atmospheric dispersion, is then necessary. Whereas carbon dioxide is often stored and handled under 2-phase or supercritical conditions associated to storage under significant pressure, CO2 ice formation is possible in case of accidental leak due to specific properties under atmospheric conditions due to its triple point. Then, this CO2 flakes creation may be followed by the liquid/solid CO2 pool formation on the ground. Due to the important pressure drop of the fluid during the leak, a significant expansion phase, a high rate of air entrainment and a huge temperature drop followed by a dense and cold cloud formation should be considered. These several steps described just above may be considered for releases of other substances stored under 2-phase conditions like ammonia and chlorine. But, due to a really high storage pressure and a specific triple point, it seems interesting to look at the models used for classical 2-phase releases to make sure that these ones remain adapted to this new situation. First, a glance at phase and enthalpy diagrams will let distinguish specific thermodynamic properties of CO2 and unsolved questions about carbon dioxide releases in the atmosphere as well. Then, the major focus of this document will consist in analyzing the different ways of determination of the liquid flash fraction χ. This fraction lets assess the rate of liquid carbon dioxide which becomes gaseous just after the point of release. This survey will specially let identify if these methods of determination (isenthalpic or isentropic methods), initially adapted to other substances, fit well to carbon dioxide. And the last but not least, we will see the influence of the χ value on the atmospheric dispersion of the CO2 cloud in terms of size and density.

PHASE DIAGRAM AND ENTHALPY DIAGRAM

Before achieving some calculations, a first insight of carbon dioxide properties looks interesting. The specific properties of carbon dioxide (CO2) are well illustrated by the two following diagrams, i e the phase diagram and the enthalpy diagram [START_REF] Ullmann | Ullmann's encyclopedia of industrial chemistry[END_REF]: These diagrams let distinguish the triple point appears at -56.4 °C and 4.1 barg. Below this point, carbon dioxide is either gaseous or solid. Due to these features, carbon dioxide can be stored and handled under high pressure (20 barg, 50 barg or more) under "dense conditions". It gathers 2 different situations: either the liquid and gaseous phases are present together in the storage or pressure and temperature in the storage are sufficiently high to make one unique supercritical phase. So, storage of carbon dioxide under dense phase conditions produces specific atmospheric releases where three phases of CO2 appear. Respectively, for 2-phase and supercritical storages:

1. the substance is initially liquid and gaseous, due to initial storage conditions; CO2 is then solid and gaseous in atmosphere, 2. the substance is initially supercritical due to initial storage conditions; CO2 is then solid and gaseous in atmosphere.

The rest of the document is focused on an academic study of atmospheric release of carbon dioxide. This case regards carbon dioxide stored under 2-phase conditions, i-e 57 barg and 293 K and released in the atmosphere. Just the liquid part of the 2-phase storage is studied. Carbon dioxide is assumed to be still liquid at the point of release and to become partly gaseous thereafter during the expansion zone. The question concerns the assessment of the rate of liquid carbon dioxide which becomes gaseous in the expansion zone. In this case, carbon dioxide is initially liquid and becomes solid and gaseous in the expansion zone. To make the document clear, liquid and solid heat capacity of CO2 are represented by Cp; in the same way, heat of vaporization and heat of sublimation of CO2 are assigned by ∆H.

Similarly, boiling point of carbon dioxide is the same parameter as sublimation point.

QUESTIONS

The important pressure drop of the fluid and the creation of solid carbon dioxide, during the leak, raise several matters which differ from classical 2-phase release like ammonia or chlorine. Especially in terms of carbon dioxide cloud formation and dispersion in atmosphere thereafter. The parameters concerned are temperature, density, volume of the carbon dioxide cloud and the air entrainment rate, just after the expansion zone, as well.

CO2 ice formation

The formation of ice, inside the pipe, appears possible if the pressure drop is sufficient to go below the pressure of 4.1 barg. Some experiments enlighten the possibility to form ice in different situations like through a safety valve on a CO2 refrigeration system [START_REF] Huang | Experimental study of throttling of carbon dioxide refrigerant to atmospheric pressure[END_REF][START_REF] Huang | Theoretical analysis of deposition and melting process during throttling high pressure CO 2 into atmosphere[END_REF][START_REF] Liu | Agglomeration process of dry ice particles produced by expanding liquid carbon dioxide[END_REF]. But at this stage, it remains really difficult to estimate the amount of ice CO2 formed for several reasons:

3. the pressure drop inside the pipe which depends significantly on the situation of leak (this problem is true as well for the 2-phase classical releases) [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]7],

4. the difficulty to form solid particles of CO2 inside a flow under high pressure [START_REF] Huang | Experimental study of throttling of carbon dioxide refrigerant to atmospheric pressure[END_REF][START_REF] Huang | Theoretical analysis of deposition and melting process during throttling high pressure CO 2 into atmosphere[END_REF],

5. the possibility to blow away layers of ice particles with high pressure [START_REF] Huang | Experimental study of throttling of carbon dioxide refrigerant to atmospheric pressure[END_REF][START_REF] Huang | Theoretical analysis of deposition and melting process during throttling high pressure CO 2 into atmosphere[END_REF],

6. the possibility to agglomerate solid particles of CO2 [START_REF] Huang | Experimental study of throttling of carbon dioxide refrigerant to atmospheric pressure[END_REF][START_REF] Huang | Theoretical analysis of deposition and melting process during throttling high pressure CO 2 into atmosphere[END_REF],

7. the interaction between the liquid phase and the solid phase, especially the possibility to melt the solid particles of CO2 with the liquid [START_REF] Huang | Experimental study of throttling of carbon dioxide refrigerant to atmospheric pressure[END_REF][START_REF] Huang | Theoretical analysis of deposition and melting process during throttling high pressure CO 2 into atmosphere[END_REF].

Liquid flash fraction in the expansion zone

The expansion zone, located just after the point of release, has been specifically studied in this document. In this zone, a part of the liquid carbon dioxide released, flashes in the atmosphere to become gaseous; and the rest becomes solid. Thermodynamic considerations have been gone through in the following paragraphs to explain the physics implied in this zone. The parameter χ is the mass rate of the liquid carbon dioxide which becomes gaseous at the end of the expansion zone in the atmosphere. This zone where flashing occurs is still studied for 2phase substance like ammonia or LPG [START_REF] Bonnet | Experimental Study of Accidental Industrial LPG Releases: Rain out investigation[END_REF]. INERIS moves on working on this topic [START_REF] Witlox | Two-phase jet releases and droplet dispersion: scaled and large-scale experiments, Droplet size correlation development and model validation[END_REF].

Airborne solid rate

Beyond this step, the ice carbon dioxide, formed during the stages before, becomes airborne or falls on the ground like "rainout". The airborne part progressively sublimates, due to its warming by the air entrained in the carbon dioxide cloud, and creates gaseous CO2. The "rainout" part either sublimates from the ground and feed the CO2 cloud or stays trapped on the ground without providing gaseous for the cloud or can be re-entrained in the jet [START_REF] Liu | Agglomeration process of dry ice particles produced by expanding liquid carbon dioxide[END_REF].

The distribution and the behavior of the different forms of solid carbon dioxide are currently difficult to evaluate. The way of release the carbon dioxide is probably crucial as well. A distinction should be made between massive instantaneous releases and continuous releases [START_REF] Witlox | Modelling of discharge and atmospheric dispersion for carbon dioxide releases[END_REF].

In this document, a specific focus has been carried out concerning the liquid flash fraction χ determination. The goal of this stage is to estimate the parts of solid and gaseous carbon dioxide at the end of the expansion zone. The formulas, adapted to 2-phase ammonia or chlorine releases, have been used to assess its relevancy in case of carbon dioxide leak.

LIQUID FLASH FRACTION χ χ χ χ

In the following steps, supercritical storage has not been specifically studied. The 2-phase storage has only been focused on.

As described before, the expansion zone regards the place where the liquid carbon dioxide, initially stored under 2-phase conditions, is released until the atmospheric pressure. In this zone, a part of the liquid released, flashes to become gaseous. Indeed, this flash is due to the pressure drop from the storage pressure to the atmospheric pressure which makes the liquid superheated. This specific state matches the situation under which liquid temperature is superior to its boiling point at 1 barg. The beginning of the expansion zone is the point of release where the CO2 is still liquid. The end of this zone matches location where the substance is under atmospheric pressure at its boiling/sublimation temperature. The rate of the liquid vaporized χ in this zone may be estimated with several methods [Britter,[START_REF] Britter | Dispersion of two phase flashing releases, Fladis field experiment[END_REF] . Some of them are described in the present document. These ways of determination are applied for the cases of classical 2-phase releases for ammonia, chlorine or LPG for instance. The storage pressure for all of these substances does not exceed 7 or 8 barg in most situations (storage, transport). These values should be compared to the carbon dioxide; for instance, its saturated vapor pressure is about 57 barg at 293 K. So, considering this last point and the specific triple point of the substance, it appears logical to wonder if the methods used to estimate liquid flash fraction χ remain true.

To compare the methods results in terms of liquid flash fraction χ, one academic case has been analyzed: the initial point is the point of release where CO2 is 100 % liquid and the final stage fits to the situation where CO2 is partly gaseous and partly solid at boiling temperature under atmospheric pressure. This case fits to either continuous releases or instantaneous releases.

Determination of χ for an isenthalpic release (Cp and ∆Η constant)

One way to determine the flash fraction χ  [START_REF] Britter | Dispersion of two phase flashing releases, Fladis field experiment[END_REF] consists in assuming an isenthalpic release from the point of release to the end of the expansion zone and to link the temperature drop of the liquid to the mass of liquid vaporized. In this case, it appears possible to consider the following equation:

χ d * H ∆ (Cp(-dT) * ) χ - (1 = (1)
χ : flash fraction (-),

Cp: heat capacity of the liquid (J/kg.K), dT: temperature drop of the liquid (K),

∆Η: liquid heat of vaporization (J/kg), dχ: mass fraction of gas vaporized (-)

In this expression, it is clear that Cp represents either the liquid capacity or the solid capacity depending on the temperature of the carbon dioxide. It is similar for ∆Η. fraction χ is determined by considering the energy amount available in the liquid phase to vaporize a part of the liquid due to the temperature difference between the storage temperature and the boiling temperature of the liquid under atmospheric pressure. Indeed, the liquid part of the 2phase substance appears in a superheated state once released in the atmosphere due to a temperature above the boiling temperature. By integrating the equation ( 1), from the storage conditions to the atmospheric conditions after the flash zone, the expression for χ can be written:

H T i ∆ / ) - Cp(T b e - 1 χ = (2)
T b : liquid temperature at the end of the expansion zone, i-e the boiling/sublimation point (K), T i : initial liquid temperature (K) considered as storage temperature.

Discussion about ∆Η and Cp variations

To enlighten the importance of assuming ∆Η and Cp constant, the carbon dioxide release has been compared to two other substances which can be stored under 2-phase conditions at ambient temperature: ammonia (NH3) and chlorine (Cl2). ∆Η and Cp have been chosen:

1. first at 293 K, the storage temperature,

Influence of ∆Η and Cp variations on the determination of χ

In order to assess influence of Cp and ∆Η on the flash fraction determination, three sets of data, (Cp1, ∆Η1), (Cp2 and ∆Η2) and (Cp3 and ∆Η3) have been used. The flash fraction calculations appear in the following tables: Table 3: final results regarding the flash fraction with the three sets of data couple

The flash fraction remains roughly constant using the two different sets of data for chlorine and ammonia releases. In the three cases, the flash fraction is estimated between 16 % and 17 % for chlorine and between 16 % and 19 % for ammonia. So, the liquid flash fraction estimation from the equation (1) seems to be adapted. Concerning the carbon dioxide, the flash fraction varies between 23 % (at 195 K) and 93 % (at 293 K), via 45 % (average). This huge range is due to strong variations of the values of the data couple (Cp, ∆Η). This statement seems to indicate the choice of the values for Cp and ∆Η at 293 K, at 195 K or averaged has a great influence of the flash fraction determination with the equation [START_REF] Harper | The inclusion of CO2 as a hazardous substance in the Seveso Directive[END_REF]. Assuming an isenthalpic release and the couple (Cp, ∆Η) constant is not efficient to precisely estimate χ  regarding CO2. To confirm this statement, another way of calculations has been selected by assuming an isentropic release with (Cp, ∆Η) constant.

Determination of χ for an isentropic release (Cp and ∆Η constant)

Assuming an isentropic release is described by [START_REF] Britter | Dispersion of two phase flashing releases, Fladis field experiment[END_REF][START_REF] Van Den Bosch | Methods for the Calculation of Physical Effects, TNO, "Yellow Book[END_REF]. The useful equation to estimate χ in this case is:

) ln( T * Cp χ b b i T T H ∆ = (3)
But once the calculations carried out, the statement is the same as before: the flash fraction varies a lot. More precisely from: 4. 21 % with Cp and ∆Η determined at 195 K, 5. more than 100 % with Cp and ∆Η determined at 293 K, 6. 49 % with averaged Cp and ∆Η from 293 K to 195 K. Similarly, considering an isentropic release with the data couple (Cp, ∆Η) constant is no more efficient due to too much influence on the flash fraction determination. Considering an isenthalpic or isentropic release by using a constant couple (Cp, ∆Η) is not appropriate. Before finding out a solution, it is interesting to focus on the enthalpy diagram to understand why common equations, used in the case of ammonia, are not adapted to carbon dioxide.

Illustration of the ∆Η and Cp variations

In order to illustrate the difference of properties concerning enthalpy of vaporization ∆Η between ammonia and carbon dioxide, a comparison of ∆Η have been gone through from these graphs [START_REF] Kirk-Othmer | Encyclopaedia of Chemical Technology[END_REF][START_REF] Ullmann | Ullmann's encyclopedia of industrial chemistry[END_REF]: For instance, from 8 bar abs. to 1 bar abs., enthalpy of vaporization is roughly constant regarding ammonia. Indeed, a variation from 1350 kJ/kg to 1200 kJ/kg seems appropriate. But unlike ammonia, from 57 bar abs. to 1 bar abs., enthalpy of vaporization varies too much to be considered constant for carbon dioxide. Besides the toughness to assess difference between ∆Η at 293 K and ∆H at 195 K, the comparison between ∆Η at 293 K (150 kJ/kg) and ∆Η at 220 K (350 kJ/kg) already indicates a significant variation for ∆Η. Similarly, the heat capacity Cp variations look significant from 293 K to 220 K as illustrated by the figure 3 B). Once ∆Η and Cp variations of CO2 illustrated, it is interesting to determine the flash fraction χ by using a different way from equations ( 1) and (3).

Another method to determinate χ : the enthalpy diagram

As enlightened before, heat capacity Cp and ∆Η vary too much to be assumed constant from 293 K to 195 K. Considering this point, the use of the enthalpy diagram has been achieved. The diagram lets estimate Cp and ∆Η variations. By assuming an isenthalpic path from the initial to the final conditions, the liquid flash fraction χ can be determined. The initial point of the thermodynamic path considered is: 1. the point of release just before the expansion zone. Carbon dioxide released is liquid at the point of release but is from 2-phase storage (293 K, 57 barg). Whereas the final point of the thermodynamic path considered is:

2. the point at the end of the expansion zone. There, carbon dioxide is under atmospheric pressure (0 barg) at sublimation temperature, i e 195 K.

From the graph, initial enthalpy Hi is estimated to 470 kJ/kg (at the point of release). By considering an isenthalpic path, final enthalpy Hf is the same as Hi (at the end of the expansion zone).

From section 2 of [START_REF] Perry | Perry's Chemical Engineer's Handbook[END_REF], ∆H is close to 570 kJ/kg at 195 K under atmospheric pressure. From the graph, H vapor = 640 kJ/kg. Consequently, H solid = 70 kJ/kg. Considering all of this information, the liquid flash fraction equals:

solid vapor solid H H H - - H χ f = (4)
By applying the equation ( 4) to this case of release, the liquid flash fraction χ = 0.7, i-e 70 %. This value should be compared to 23 %, 45 % and 93 % determined by considering ∆Η and Cp constant.

DISCUSSION

As seen before, the liquid flash fraction χ values vary with the way of determination. For instance, by considering an isenthalpic release from an initial storage at 293 K, the χ value is set from 23 % to 93 % via 45 % by considering Cp and ∆Η constant at different temperatures and 70 % by using enthalpy graph (the more relevant way). Due to carbon dioxide is either solid or gaseous under atmospheric conditions, the rest of the percentage is considered solid. The solid carbon dioxide appears under small ice particles which are either airborne and stay in the cloud or "rainout" and fall on the ground. At this step, it is still difficult to precisely assess the phase distribution in the cloud in terms of gaseous, solid and airborne solid carbon dioxide. Data from tests are still missing.

Continuous releases

The first way is not to consider any "rainout" due to carbon dioxide particles sublimate very quickly. This assumption is assumed to be true by [START_REF] Witlox | Modelling of discharge and atmospheric dispersion for carbon dioxide releases[END_REF] for most not very large and continuous releases. However, rainout could be more likely for large diameter and/or low pressure releases. Indeed, in this case, heat brought by the air entrained could not be sufficient to let a very quick sublimation. Moreover, tests carried out by [START_REF] Liu | Agglomeration process of dry ice particles produced by expanding liquid carbon dioxide[END_REF] for small diameter (few millimeters) releases seem to indicate that carbon dioxide ice trapped on the ground is quasi immediately re-entrained in the CO 2 jet. But in this case, the height of release is significantly weak (few millimeters). It appears necessary to carry out some experiments to go further regarding this step. INERIS has planned to achieve such a work in 2011.

Instantaneous releases

By using c and considering no "rainout", the density and the volume of the initial cloud by expansion, after flashing and before mixture with air, can be assessed by the following formulas:

sol gas cl ρ ) χ - 1 ( ρ χ ρ 1 + = (5)
ρ cl : global cloud density at 195 K (kg/m 3 ), ρ gas: gaseous cloud density at 195 K (kg/m 3 ), ρ sol : solid cloud density at 195 K (kg/m 3 ), and

cl cl cl M V ρ = (6)
V cl : global cloud density at 195 K (m 3 ), and Μ cl : global mass cloud at 195 K (kg).

The temperature of this initial cloud is assumed to be the same as the temperature at the end of the expansion zone, i e 195 K. By considering gaseous density close to 2.8 kg/m 3 , solid density about 1550 kg/m 3 and a mass released about 1000 kg (for instance), four situations should be distinguished upon the value of χ: Thereafter, the cloud features are largely modified due to its interaction with the air around. At this stage, a crucial point, which needs to be focused on in a future work, is enlightened: the air entrainment rate following the expansion zone. This crucial step constitutes the link between the source term, discussed in this document, and the last step regarding the atmospheric dispersion of the CO2 cloud.

The estimation of the air entrainment rate, just after the expansion zone, is especially interesting for two main reasons:

the high concentrations of interest about 100 000 ppm. Such a level leads to focus on a field located really close to the point of release where the air entrainment may have a strong influence in terms of concentration value, the calculations of carbon dioxide cloud enthalpy which influences the rates of airborne solid and gaseous carbon dioxide; with a consequence on cloud temperature, density and volume. These parameters are necessary to estimate in order to carry out the atmospheric dispersion of the cloud, denser than the air, which moves along the ground with accumulation of the substance [START_REF] Riou | Comment estimer la dispersion des gaz lourds, EDF[END_REF] and [START_REF] Markiewicz | Mathematical modelling of the heavy gas dispersion, the faculty of environmental engineering[END_REF].

CONCLUSION

In the future, significant industrial facilities could handle and transport important amount of carbon dioxide due to global warming considerations. The CO2 could be used under different conditions and in particular under 2-phase conditions. In this situation, an accidental leak of carbon dioxide would provide an expansion zone close to the point release due to significant depressurization from the storage conditions to atmosphere. In this specific zone, a part of the liquid carbon dioxide flashes to become gaseous. And due to carbon dioxide features under atmospheric conditions, the rest of the liquid becomes solid. From this study, it would appear that common formulas used to assess the rate of the liquid flashed in the expansion zone, for ammonia or chlorine for instance, are not so appropriate. This statement is due to several topics. First, CO2 depressurization, in case of leak, is really more significant than in other cases like ammonia or chlorine. Secondly, the vaporization enthalpy ∆Η and the liquid heat capacity Cp vary too much during the depressurization to be considered constant. And finally, the triple point of CO2 at 4.1 barg and 195 K creates discontinuity for ∆Η and Cp with the formation of solid and gas from a liquid release. So, to get over this matter, it appears interesting to consider the enthalpy diagram to estimate precisely the liquid flash fraction. This fraction has an influence on the formation of the initial cloud just after the expansion zone in terms of density and volume. Thereafter, the crucial point which links the initial cloud after expansion and the atmospheric dispersion is the air entrainment phase.

Figure 1 :

 1 Figure 1: phase diagram (A) and Enthalpy-Pressure-Temperature-Entropy diagram (B) of the carbon dioxide from the Ullmann's encyclopedia vol.A5, 167

Figure 3 :

 3 Figure 3: Mollier diagram for ammonia (A) and for carbon dioxide (B)

Figure 4 :

 4 Figure 4: thermodynamics way of the 2-phase carbon dioxide from the storage conditions to the end of the expansion phase

Table 1 :

 1 2. then at the boiling or sublimation temperature of each substance. The boiling points of chlorine and ammonia at ambient pressure are 239 K, 240 K whereas the sublimation point of carbon dioxide at ambient pressure is 195 K. Concerning the carbon dioxide, ∆Η and Cp have been determined at 195 K by taking into account the following graphs. These graphs let assess the liquid phase of the carbon dioxide from 293 K to 216 K and the solid phase of the carbon dioxide from 216 K to 195 K:3. and with averaged ∆Η and Cp between 290 K and the temperature of vaporization/sublimation. Averaged values for  ∆Η and Cp have been set up by considering the graphs above as well. The following table regards Cp and ∆Η of the three substances:

		Cp1 at 293 K (J/kg/K)	Cp2 at Tb/sub (J/kg/K)	Averaged Cp3 (J/kg/K)	Cp1/Cp2 (-)	Cp1/Cp3 (-)	293 K ∆Η1 (J/kg) at	Tb/sub ∆Η2 (J/kg) at	(J/kg) Averaged  ∆Η3 ∆Η1/∆Η2 ∆Η1/∆Η3
	NH3	4660	4480		4560	1,04		1.02	1.19E06 1.37E06 1.28E06	0,87	0,93
	Cl2	900	930		915	0,97		0.98	2.55E05 2.87E05 2.70E05	0,89	0,94
	CO2	4100	1500		2150	2,73		1,90	1.53E05 3.71E05 3.50E05	0,41	0,45
					Cp 1 /∆Η 1 (1) Cp 2 /∆Η 2 (2) Cp 3 /∆Η 3 (3) (1)/(2)	(1)/(3)
			NH 3	3.92E-03	3.27E-03	3.56E-03	1,20	1.10
				Cl 2	3.52E-03	3.23E-03	3.39E-03	1,09	1.04
			CO 2	2.68E-02	4.04E-03	6.71E-03	6.63	4.34

Cp and ∆Η variations

Whereas the ammonia and the chlorine liquid heat capacity remain roughly constant from the storage temperature to the atmospheric temperature, the solid/liquid heat capacity of carbon dioxide decreases significantly with temperature decrease. In the same way, the heat of the liquid sublimation/vaporization of CO2 varies a lot with temperature in comparison of ammonia or chlorine. To go further, the ratio (Cp / ∆Η), used to assess χ, is compared as well. Results are summarized in the following table:

Table 2 :

 2 ratio (Cp / ∆Η) variations with temperature variations Figure2: heat of sublimation ∆Η[START_REF] Perry | Perry's Chemical Engineer's Handbook[END_REF] and liquid heat capacity Cp[START_REF] Witlox | Modelling of discharge and atmospheric dispersion for carbon dioxide releases[END_REF] from 290 K to 190 K As expected, values reported in table 2 further reveal that temperature slightly affects ammonia and chlorine properties while, in the case of carbon dioxide, the ratio Cp/∆Η increases more than six-fold as temperature changes form atmospheric to boiling/sublimation value. On the whole, liquid/solid heat capacity Cp and heat of the vaporization/sublimation ∆Η of the carbon dioxide cannot be considered constant from 293 K to 195 K. Indeed, Cp value at 195 K is roughly twice inferior to the value at 293 K. The ∆Η value at 195 K is more than twice superior to the value at 293 K. In the same way, the ratio (Cp / ∆Η) value at 195 K is nearly seven times inferior to the value at 293 K. And the last but not least, the choice of averaged values, between 293 K and the boiling or sublimation temperature, for Cp and ∆Η, provides different value regarding the ratio. Consequently, liquid/solid heat capacity Cp and heat of vaporization/sublimation ∆Η should not be considered constant for the carbon dioxide during the release. So, the next step consists in assessing the influence of the Cp and ∆Η variations on the liquid flash fraction determination χ.

Table 4 :

 4 density and volume initial cloud for 1000 kg released Values in the table do not take into account any solid CO2 rainout, However, in case of massive instantaneous release; it is possible to imagine a significant part of the carbon dioxide ice can be trapped on the ground. Experimental data is missing at this stage to confirm this statement. Moreover, density and volume of the initial cloud are just relevant after the expansion zone.

	Gaseous CO 2 (χ)	Solid CO 2 (Airborne+rainout)	Density of the initial cloud (kg/m 3 )	Volume of the initial cloud (m 3 )
	23 %	77 %	12.1	80
	45 %	55 %	6.2	160
	93 %	7 %	3	330
	70 %	30 %	4	250


Frédéric Tognet

Cyrille Turmeau

Ti L Ha

Eric Tarnaud

Laurence Rouïl

Bertrand Bessagnet

Enric Robine

Yannick Morel

NUMERICAL MODELLING OF MICROORGANISMS DISPERSION IN URBAN AREA: APPLICATION TO LEGIONELLA

Keywords: Legionella, microorganism, biological aerosols, dispersion, numerical modelling, biological model, biosamplers

Dispersion modelling is often used to estimate potentially contaminated areas in case of accidental release of microorganisms in the atmosphere. In the specific case of Legionella, accidental spread in the atmosphere due to contaminated cooling towers system may occur over distance larger than 10km. In addition, most cooling towers are located in urban areas where dispersion due to obstacles is complex. In this case, dispersion models have to take into account complex flows and microphysical processes that occur within the plume and may have an impact on the survival of the microorganisms. To estimate the concentration of microorganisms in these areas, a specific module has been developed within the lagrangian dispersion model Micro Swift Spray (MSS, Aria technologies). This module takes into account microorganisms outside or inside water liquid droplets and microphysical interaction inside the plume. A simple biological module governing the survival of airborne microorganisms has also been implemented in the dispersion model. In order to evaluate this model, a field campaign of biological aerosols dispersion was performed by CSTB (Champs-sur-Marne, France) on June 23 rd , 2009. Spores of bacillus atrophaeus (usually referred to as BG) initially contained in a water tank were disseminated in a suburban area from a source at 3 meters above ground level. Air was sampled by DGA MNRBC (Vert-Le-Petit, France) at 4 various locations from 50 to 300 meters from the source to monitor NG concentration. Direct impaction onto Petri dishes was performed with slit-samplers and sixstage Andersen impactors. Wetted-wall cyclones and SKC Biosamplers were also implemented in order to sample air and generate liquid samples. Dispersion modelling for this campaign has been carried out by INERIS using the microorganism module developed in MSS. The results show that predicted concentrations and in situ measurements are in agreement. MSS Model was implemented to simulate legionella airborne dispersion from a virtual cooling tower at the same location. The biological model has been activated. Results show that the impact of biological model on airborne concentration is significant.

INTRODUCTION

Legionella is an opportunistic bacterium responsible for Legionnaires' disease, a respiratory infection for humans. The occurrence of Legionella community has sometimes been linked with transmission of an infectious agent from cooling towers (CT). During the episode that occured in Pas-de-Calais (France) in 2003, a transport of Legionella up to a distance of 12 km has been observed [START_REF] Nguyen | A community-wide outbreak of legionnaires disease linked to industrial cooling towers--how far can contaminated aerosols spread[END_REF]. In such an incident, the use of a numerical dispersion model would be helpful to test the influence of several potential sources to guide the search of the CT involved as well as to estimate the extent of the potentially contaminated area.

Although the description of the CT plumes is conventionally based on the use of numerical models; their effectiveness to assess or anticipate the impact of dispersion on the local population's exposure remains to be validated. Moreover, just a few models allow the simulation of the dispersion of microorganisms. For instance, [START_REF] Lighthart | Simulation of airborne microbial droplet transport[END_REF], introduced a bioaerosol droplet dispersion model with evaporative processes and biological decay in the aerosol droplets, to obtain a reasonable fit to observed data from a field spray of a General Ecosystem Model within the first 30 meters from the source [START_REF] Ganio | A comparison between computer modeled bioaerosol dispersion and a bioaerosol field spray event[END_REF]. In the specific case of Legionella, accidental release in the atmosphere due to contaminated cooling towers may occur over large distance. Moreover, most CTs are located in urban or sub-urban areas where accounting for the dispersion due to obstacles is challenging. In this case, a dispersion model has to take into account complex flows and microphysical processes that act within the plume. That is why a specific module was developed inside the lagrangian model MSS [START_REF] Tinarelli | Micro-swift-spray (mss): a new modelling system for the simulation of dispersion at microscale. General description and validation[END_REF]. The objective of this study supported by the AFSSET( Agence Française de Sécurité Sanitaire de l'Environnment et du Travail) consisted in estimating the ability of this tool to predict the potential spread of microorganisms by an anthropogenic source resembling a CT, by comparing data from the simulation with experimental measurements. This work refers only to the transport of bioaerosols in the air. The numerical model and the modelling methodology is first presented, then the experiments, which consisted of in situ dispersions by a biological tracer (spores of Bacillus atrophaeus) on a suburban site. Results and comparison are presented. An approach of the dispersion modelling of legionella airborn from a virtual CT is made. The work was broken down as follow: the DGA has provided measurements of concentration and bioaerosol particle size which were then compared with the results of simulations performed by INERIS. CSTB provided the technology for assembling the source, some measurements concentration as well as the set up of the biological model.

NUMERICAL MODEL

In this study, the lagrangian model MSS (Aria technologies) was chosen. The wind field is reconstructed from average measured variables (temperature, relative humidity, wind direction and intensity) with the meteorological preprocessor (Micro Swift). The final wind field (which is mass consistent) is diagnosed using analytical corrections induced by the presence of obstacles, considering the topography and buildings. Thus, the recent module "obstacle" of MSS improves the accuracy of the results for near field and urban areas. The source code must also be accessible in order to implement the modules described above which are specific to bioaerosol (cultivability of bacteria) and their sources (evaporation-condensation in the plume of a CT). Thus, a specific module was developed inside the lagrangian model MSS. These developments were conducted in partnership with ARIA Technologies. The underlying principle is that the 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 computed for each droplet and the evolution of the temperature of the droplet is driven by the generic Lagrangian equation (RS [START_REF] Miller | Evaluation of equilibrium and non-equilibrium evaporation models for many-droplet gas-liquid flow simulations[END_REF]. Water vapour diffusion into the air is solved on a eulerian grid. Finally, the evolution of the density and temperature within the plume is evaluated using the Glandening (1984) equations. 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 approach specific to L. pneumophila has been applied by CSTB. This approach is used to set a module of sustainability that is superimposed on the transport calculations.This biological model based on an exponential survival of microorganisms is applied at the specific moment when microorganisms are transported as free airborne and are directly exposed to radiation, air temperature and humidity.

EXPERIMENTAL DEVICE

Dispersion system

The dispersion system has been designed by CSTB to mimic a small CT. The emission is performed at 3 m height above the ground, using 8 high pressure nozzles located on the section of a vertical pipe mounted on a mobile scaffolding (Figure 1). The nozzles produce droplets that have a median volume diameter of 30 microns, similar to that of ART available. The device is operated for up to 30 minutes above a tank containing the spore in suspension. Figure 2. Spores of B. atrophaeus (MEB)). It is an environmental (telluric) microorganism devoid of pathogenicity or toxicity, used as environmental tracer. Another advantage of such a tracer is our ability to detect these harmless particles at very low concentrations. In confined spaces, the NMAD of one spore is 1.4 microns (measured by Andersen). Assuming the natural tendency of these spores not to aggregate in diluted suspensions, a concentration below 3,5×10 7 spores / ml lead to an average of 1 spore for 2 drops of 30 microns diameter. This induces a low probability for aggregates formation.

Biological agent

Biocollectors

Instrumentation consisted in Agar impactor (yielding particle/L air results; one such particule can contain one or more bacterial cells) or in a liquid medium (cfu/L air ; cfu : colonny-forming unit ie. viable isolated bacterial cell), arranged on a support of 80 cm. Reproducibility of these instruments in controlled conditions is usually good (relative standard deviation of the air concentration of 5% for agar impactors and 10% in a liquid medium).

 impactor slot STA 203, New Brunswick (NBSS; Fig. 3) that operates at a aspiration flow rate of 50 L / min, with a D ae 50 smaller than 0.6 microns.  Andersen 6-stage impactor TE-10-800 Tisch Environmental (AND, Fig. 4) that operates at an aspiration flow rate of 28.3 L / min with an omni-directional sampling head. There are six Petri dishes containing suitable growth medium, kept under sieves of different pore size. Each sieve has 400 pores. The D ae 50 falls from 7 to 0.65 microns on the different stages.  Soprano wet wall cyclone Bertin Technologies (CYC; Fig. 5) set at a rate of 550 L / min (± 5%). The air inlet of the sampling head is multidirectional.  Impinger Biosampler of SKC (SKC, Fig. 6): aspiration from 12.5 L / min for an initial 15 mL of liquid.

Campaign set up

A first campaign was carried out to set up the dispersion system and to check the sensitivity of the biocollectors. For this test, the source was installed on a roof. All instruments downwind and within 300 m (on roofs and on the ground) detected bioaerosols. The second campaign ran a few months later on the same day time to take benefit of similar weather conditions for the three tests (dispersions C, D and E). For the final test (dispersion E), the concentration of the spore suspension was multiplied by a factor 7. The instruments were located 20 to 300 m away from the source according to the distribution of Table 1. Points 2 and 4 were located on the roofs of the buildings.

NBss X X X X X X AND X X X SKC X X X X CYC X X X
Meteorology (3 points) a) 30m above ground and 30m away from the sourceb) 2m above ground at source locationc) 2m above ground at point 3

The final disposition was adjusted at the beginning of the testing phase (Figure 7), according to local meteorological measurements (the generator, central site and 30 m above ground).

RESULTS

Table 2 summarizes the operating conditions for the three tests on Tuesday, June 23rd, 2009, a fair day with dry air, and unstable atmosphere (Pasquill class: A). For the two concentrations of the tracer suspension and the three distances to the source, cultivable particles diameter are between 1.1 and 2.1 microns. The stability of NMAD with the distance shows that the aerosol is rapidly "dry" 20 m away from the source. Furthermore, the value of NMAD does not change when the concentration of the source is multiplied by a factor 7 (still below the theoretical limit concentration to avoid Bg aggregates). This implies an absence of aggregates. Therefore, the high value of NMAD (slightly above one spore size) could be explained by an agglomeration of inert particles on the spores in the dispersion (an effect of scavenging when the droplets are not totally dry). Otherwise, we have considered only particles containing a single cultivable spore, offering thus a good consistency between the biocollectors units (ie. 1 part./Lair = 1 cfu/Lair).

To check the model ability to predict the aerosol dilution according to the distance to the source, we compared the evolution of average concentrations from the NBSS and MSS simulations. Normalization to 10 cfu/m 3 of the point 6 results (exposed to lower turbulence) overcomes a possible systematic bias between biocollectors and simulations. For each experiment, results show a decrease in average concentrations of cultivable spores away from the source (examples in Figure 8), with a good agreement between model results and observations (discrepancies are discussed below). However, a significant difference appears at point 8 (300 meters), for which a relatively high concentration measured is not captured by the model. An accumulation due to a recirculation zone behind a building may have broken away from the obstacle and then impacted point 8. However, such a phenomenon cannot be represented by our dispersion model performed on a 1 minute stationary situations basis. To represent the dispersion dynamics, we have normalized the maximum concentrations of NBSS and the MSS model and compared their time evolutions. Only items 4 and 6 (supposed to be better captured by the model) were studied. The general dynamical model is close to experimental results (examples in Figure 9). The duration of the bacterial cloud and amplitude fluctuations are thus highly correlated. In some situations, hot concentrations may be missing or added according to field measurements. The use of stationary weather conditions averaged over one minute can be responsible for this feature. For the non-normalized values, significant differences were found between the concentrations measured by some biocollectors and simulations. The average ratios of concentrations (model / biocollector) are respectively 9 -1.9 -1.6 and 5 for the NBSS, the AND, the SKC and CYC. If the model can be questioned, an experimental bias must also be considered because of the existence of a ratio of 1 to 6 between each type of biocollector. In addition, obtaining ratios below a factor of 1.3 between the NBSS and the AND when used in a controlled environment indicates that the problem is specific to an external metrology. If the relationship between our four types of biocollectors and with the MSS model are consistent on average, the dispersion of ratios remained significantly higher than in a controlled way of testing (laboratory facility). Thus, the relationship model {/} biocollectors on these tests span an order of magnitude and differences in relative efficiencies of biocollectors reach a factor of 5 (less than 1.5 in laboratory conditions). The origins of this variability would be distributed between the model and therefore the fact of using instrumentation outdoors.

COOLING TOWER ("CT") DISPERSION MODELLING

A characterization of the source term "CT" was conducted by CSTB. It focuses on the size of the droplets and transfer of microorganisms from the water to aerosols during operation. A functional type CT was made available by the company Balticare and installed on the site of the CSTB. Its operating characteristics were controlled by prior device instrumentation. This gave information on emission parameters required for simulations. In addition, based on the work of the CSTB for the survival of L. pneumophila in aerosol according to the ambient humidity (TL Ha, 2005), two coefficients α were defined: one describing the loss of cultivability due to stress at the time of the emission (α = 0.006), the other describing the loss of cultivability during transport (α = 0.002). From the physical parameters of the CT, and the meteorological data measured during the dispersion E, we simulated a dispersion from the CT placing it at the source studied previously and using the biological module which describes the conditions of survival and cultivability of Legionella. Figure 10 shows a simulated contamination of 10 5 CFU / L within a small CT (French regulated maximum concentration value in the liquid water of a CT before its activity is stopped) where atmospheric concentrations are highest in the first 400 meters downwind of the tower.

Results show that applying the biological model leads to a cultivability of 40% at 400 meters from the source. The resulting concentration values are very low and show that contamination by isolated bacteria seems unlikely given the characteristics of CT and investigated the initial concentration used.

CONCLUSION

Qualitatively, experimental measurements and numerical results from software simulations are similar as to spatial and temporal variability of air concentrations. However, significant differences remain over short time scales. Using an Eulerian model for the diagnosis of an unsteady wind field would have limited this behavior. This remark also applies to the quantitative analysis of the results since the model predicts average concentrations about twice as high than the experimental results. Nevertheless, additional tests and sampling farther away from the source would be helpful to validate its operational use.

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 module. The type of model (Gaussian, Lagrangian, Eulerian, with specific modules or not) will depend on the expected accuracy and duration of the simulations requested by the authorities for the decision making process. Experimentally, 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.

  Spores of B. atrophaeus (known as B. globigii or BG strain reference: ATCC 9372, CIP7718, DSM 675, FRBD Batr001) are similar in size to Legionella: their average dimensions is 0.7 x 1.4 microns (Figure 1. Dispersion system.

Figure 1 .

 1 Figure 1. Dispersion system.Figure 2. Spores of B. atrophaeus (MEB)

Figure 2 .

 2 Figure 1. Dispersion system.Figure 2. Spores of B. atrophaeus (MEB)

FigureMeteorology

  Figure 3. Impacter STA203

Figure 7 .

 7 Figure 7. Modelling of the site, location of the generator and locations of collection points.

Figure 8 .

 8 Figure 8. Evolution of normalized concentrations depending on the distance to the source (disp C and E; NBSS and model)

Figure 9 .

 9 Figure 9. Time evolution of normalized concentrations (disp E--pt4 pt6; NBSS and model)

Figure 10 .

 10 Figure 10. Maximum concentration of legionella depending from the distance to the source

Table 1

 1

	: instruments location
	Collection points and distance from source	20m pt1 100m pt2 100m pt3 100m pt4 200m pt5 200m pt6 200m pt7	300m pt8
	Biosensors

Table 2 :

 2 summarize of the three dispersion experiment.

		Duration	Concentration of	dispersed	T°	RH	Pressure	Wind	Wind
		disp.	suspension	Volume				direction	intensity
	Disp C	22,5 min	4.6 x10 6 cfu/mL	6,4 L	21 °C	35 %	1019 hPa	59 deg.	2,5 m/s
	Disp D	25 min	4.3 x10 6 cfu/mL	7 L	23 °C	31 %	1018 hPa	76 deg.	3,4 m/s


B Hor

M Caudron

& M Al Heib

Experimental analysis of the impact of ground movements on surface structure

Man's activity frequently causes ground movements: the most notable examples of ground movements are provided by the mining industry in the form of subsidence. Consequently, mining subsidence can damage the surface structures. Many empirical approaches have been developed to predict structure damage using ground movements measured in greenfield conditions by assuming that the presence of a building has no effect on the settlement prediction. To improve our understanding of that effect and investigate the potential damage to the building, we have undertaken an experimental research program by means of a large three-dimensional physical model designed to simulate mininginduced ground movements. The results pointed out the effect of building structure on greenfield movements to be taken into account and enabled an analysis of the damage to the surface building, based on the existing damage estimation approaches. RÉSUMÉ L'activité humaine cause fréquemment des mouvements de terrain. Les mouvements de terrain sous la forme d'affaissement sont induits par l'exploitation minière. Les affaissements miniers peuvent endommager les structures en surface. Des approches empiriques ont été développées pour prévoir les dommages au bâti en utilisant les caractéristiques des mouvements de terrain sans prendre en compte l'effet de la présence des constructions en surface. Afin d'étudier ces effets et d'investiguer les dommages potentiels sur le bâti, nous avons réalisé une étude expérimentale à l'aide d'un modèle réduit tridimensionnel de grandes dimensions permettant de simuler des mouvements de terrain dus aux mines. Les résultats expérimentaux montrent l'influence de la présence du bâti sur les mouvements du sol, qui doit être prise en compte dans l'estimation de dommage, et permettent d'apprécier la prévision des dégradations de la structure bâtie basée sur les méthodes existantes. 1

INTRODUCTION

Subsidence of the ground surface can be regarded as ground movement which takes place due to the extraction of mineral resources. It is an inevitable consequence of mining activities and reflects the movements which occur in the mined out area. Unfortunately, subsidence can and does have serious effects on surface structures and services. Prior to construction, the ground movements should be predicted and their consequences on the structures and services should be assessed. Many research projects have been focused on the prediction of structure damage prediction by applying empirical (NCB 1975;Wagner & Schümann 1991) or analytical approaches [START_REF] Burland | Settlement of buildings and associated damage[END_REF][START_REF] Boscardin | Building response to excavation-induced settlement[END_REF][START_REF] Burland | Assessment of risk of damage to buildings due to tunneling and excavations[END_REF]). Yet, these approaches present an important limitation, since they are usually based on greenfield predictions without taking the influence of building stiffness or soil-structure interaction into account. Research using numerical analysis has allowed the effect of building stiffness on greenfield movements to be taken into account (Potts & Addenbrooke 1997). However, there is still great uncertainty about how to estimate the building stiffness itself and how to assess the soil-structure effects. [START_REF] Standing | Impact of underground works on existing infrastructure[END_REF] and [START_REF] Dimmock | Effect of building stiffness on tunneling-induced ground movement[END_REF] have provided some evidence of the effect of building stiffness and soil-structure interaction on settlement prediction for tunnelling based on observation data. [START_REF] Standing | Impact of underground works on existing infrastructure[END_REF], on the other hand, suggested that the methodology for categorising and assessing damage (according to [START_REF] Burland | Assessment of risk of damage to buildings due to tunneling and excavations[END_REF] and incorporating building stiffness into the prediction of building deformations (according to Potts & Addenbrooke 1997) developed for tunneling-induced ground movements could be applied for cases involving mining-induced subsidence if the greenfield movements can be predicted reliably.

The aim of our research is to use an experimental modeling approach by means of a large 3D physical model developed by "Institut National de l'Environnement industriel et des Risques (INERIS)". The model makes it possible to represent a variety of ground movement geometries, and then to study the impact of ground subsidence due to mining on surface buildings by taking soil-structure interaction into account. Another objective is to evaluate the building damage based on the existing approaches.

This paper first briefly outlines the characteristics of mining-induced subsidence and its consequences on the surface structure, and then summarises the damage estimation methods for buildings located in the mining subsidence area. An important development to predict ground movements by taking the stiffness of a building into account (Potts and Addenbrooke, 1997) will be also briefly presented. We will then give a presentation of the physical model and the building model representing an individual house. Much of the content of this paper is drawn from the discussion on the effect of soil-structure interaction and the estimation of building damage using data from experimental tests.

MINING SUBSIDENCE AND STRUCTURE DAMAGE

The subsidence parameters from underground excavation that are used to assess the impacts of subsidence on building structures are shown in Figure 1. When the ground movement occurs, the overlying buildings or structures are in general affected by the vertical and horizontal displacements.

Figure 1. Subsidence parameters and structure's behaviour for different positions at ground level Vertical displacement (S v ) due to mining generally results in a "subsidence trough." The maximum subsidence (S v,max ) is typically found around the trough center and the extent of the trough is limited by the angle of draw. The impact of subsidence on the building results from the occurrence of differential subsidence. The differential vertical movements give rise to ground slope (α).

Horizontal displacements (S h ) due to mining subsidence occur in such a way that points on the surface generally move towards the centre of the subsidence trough. Differential horizontal displacements give rise to horizontal strains (ε h ).

Ground strain is one of the major causes of impact to structures. Another major cause is curvature resulting from differential slope. Within the subsidence trough, convex or hogging curvature is accompanied by tensile strain (ε ht ) and concave or sagging curvature is accompanied by compressive strain (ε hc ). The hogging zone is separated from the sagging zone by the inflexion point where the maximum ground slope (α max ) and maximum horizontal displacement (S h,max )-take place. Both tensile and compressive strains can cause cracking in a building structure, but tensile strains are more difficult to accommodate since almost all components of a structure are weaker in tension than in compression.

The transfer of ground strains into the structure occurs through friction on the underside of the foundations and ground pressure on the sides of the foundations or on buried walls. The transfer is thus dependent upon the configuration and type of foundation and its orientation to the subsidence trough. The transfer of strain is also dependent upon the types of soil that are immediately below the foundation.

BUILDING DAMAGE ESTIMATION

Clearly, if an estimation of the building damage due to ground movement is to be made, the classification of damage is a key issue.

Classification of damage

Many studies have been largely carried out in the UK to classify the building damage subjected to tunnellinginduced ground movements (according to [START_REF] Standing | Impact of underground works on existing infrastructure[END_REF]. Three broad categories of building damage can be defined: (a) visual appearance or aesthetics, (b) serviceability or function, and (c) stability. As foundation movements increase, damage to a building will progress successively from (a) to (c). [START_REF] Burland | Settlement of buildings and associated damage[END_REF] divided these broad categories in six categories of damages, numbered from 0 to 5 in increasing severity (negligible to very severe damage). Normally categories 0, 1 and 2 relate to aesthetic damage, 3 and 4 relate to serviceability damage and 5 represents damage affecting stability.

Damage estimation methods

The damage criterion used in the above-mentioned classification system is based on the ease of repair of the visible damage (e.g. crack width of the brick or masonry structure). This damage criterion could be expressed in terms of the limiting tensile strain (ε lim ) developed by [START_REF] Burland | Settlement of buildings and associated damage[END_REF]. [START_REF] Boscardin | Building response to excavation-induced settlement[END_REF] went on to analyze this concept and finally showed that the categories of damage could be broadly related to ranges of ε lim as presented in Table 1.

The most used method established by the UK National Coal Board (NCB 1975) makes it possible to estimate building damage as a function of the structure length and the horizontal ground strain. This method assumes that the building has no effect on the settlement prediction and the damage parameters are calculated using the predicted greenfield movements. The NCB's approach has been proposed for very long structures. In the case of individual house considered in the present analysis, the NCB's approach is useless, since it is too imprecise. Therefore, [START_REF] Wagner | Surface effects of total coal-seam extraction by Underground mining methods[END_REF] developed another chart using the same damage parameters and based on structure damage observed mostly in South Africa. This chart, presented in Figure 2, conforms to NCB's for lengths of structure longer than 45 m.

As mentioned above, following the concept of limiting tensile strain established by [START_REF] Burland | Settlement of buildings and associated damage[END_REF], [START_REF] Burland | Assessment of risk of damage to buildings due to tunneling and excavations[END_REF] showed how, by combining the deflection ratio (Δ/L where Δ is the deflection over the structure length L) with the horizontal strain simple interaction diagrams relating (Δ/L)/ε lim to ε h /ε lim for various values of L/H could be developed. Figure 3 shows one such diagram for the case of an isotropic beam with L/H = 1 undergoing hogging. The limiting tensile strains have been converted into the associated categories of damage given by Table 1. to ε lim in the range 0.15 -0.3% as moderate to severe. However, none of the cases quoted by them exhibit severe damage for this range of strains. There is therefore no evidence to suggest that tensile strains up to 0.3% will result in severe damage.

Influence of building stiffness on settlement prediction

The prediction of ground movements due to tunnelling and mining is normally done using empirical relationships which are based on measurements from greenfield conditions. The presence of a building is assumed to have no effect on the settlement prediction and any damage parameters are calculated using the predicted greenfield movements. This is clearly an over-simplification of reality as the stiffness of the building affects its deformation.

The Potts and Addenbrooke (1997) method addresses the influence of a surface structure on the ground movements due to tunnelling and has proved to be a useful tool in the prediction of potential building damage. A series of 2D numerical analyses of tunnel construction in greenfield conditions and beneath surface beams of varying linear elastic stiffness, oriented transverse to the tunnel, were performed to assess the difference between greenfield ground movements and those modified by the structure. The building deformation parameters of the deflection ratio and horizontal strain were expressed as a fraction of those obtained for greenfield conditions. The modification factors are defined separately for deflection ratios in sagging and hogging (see Figure 4), and for horizontal strain in compression or tension, and are as follows: A design approach for predicting building deformation was proposed, which is based on factors of relative structure-soil stiffness. The relative bending stiffness, ρ* and relative axial stiffness, α*, are defined as follows:

B E 2EA α* ; B E 16EI ρ* s 4 s [5,6] Building length (m)
where B is the building width, E s is a representative soil stiffness, and E, A and I are an equivalent Young's modulus, cross-sectional area and moment of inertia for the beam. It should be noted that ρ* is not dimensionless, but rather has the dimensions m -1 in plane strain; this is discussed by [START_REF] Franzius | The response of surface structures to tunnel construction[END_REF] who proposed a revised dimensionless form for ρ*.

Based on the numerical analyses for 5-storey, 3-storey and 1-storey buildings and for a single slab, Potts and Addenbrooke presented a series of upper bound design curves for the variation of modification factor with relative bending stiffness and relative axial stiffness, each curve representing a different eccentricity ratio (e/B) (see Figure 5). A modification factor of 1 or greater represents fully flexible behaviour, because the building follows the greenfield settlement profile, whereas a modification factor close to zero represents very stiff behaviour (the deflection ratio of the building being almost zero). This paper intends to shed some light on the soilstructure interaction due to mining subsidence by means of physical modeling. To study that phenomenon we need to have a physical model that is capable of simulating the ground movements, a soil model that is a natural soil, and a building model.

The 3D physical modelsimulator of ground movements

The details of the hypotheses and conception, the constraints and limitations, as well as the monitoring method of the small-scale physical model can be found in [START_REF] Caudron | A large 3D physical model: a tool to investigate the consequences of ground movements on the surface structure, 14[END_REF]. The main hypothesis of the model is that it is designed to be used in a 1g (earth gravity) environment and not in a centrifuge. It is therefore difficult for quantitative interpretations of the experimental data to be made. The geometry scale factor may range up to 1/50. With a size of 2 m by 3 m in the horizontal plane for a height up to 1 m, we would be able to represent a soil block as large as 100 m by 150 m with a height of 50 m (Figure 6).

Considering underground mines as cavities to be the origin of the ground movements caused by their collapse, the maximum depth of 50 m may appear to be an important limit to this model. However, we chose not to model the cavity itself, but only to create the subsidence trough at ground surface, which is similar to subsidence induced by mining. This is achieved by vertically moving an "electric jack" placed at the bottom (initial state) of the model downwards to the desired opening depth (final state). This mechanism for creating the subsidence trough at ground surface is presented in Figure 7.

In our study, a geometry scale factor of 1/40 with an overburden of 0.3 m and a limit subsidence extent of 1.3 m was chosen, as it provides a geometry equivalent to a real mining case. The cross-section in horizontal plane of the jack is limited to 0.25 x 0.25 m², corresponding to 10 x 10 m² at prototype scale. The selected Fontainebleau sand is very fine and clear. The diameter of the grains varies from 0.1 to 0.3 mm with D 50 approximately 0.2 mm. A study of its mechanical properties was carried out prior to the experiments. The estimated properties of the soil mass model are presented in the Figure 7.

A measurement technique known as Digital Image Correlation (DIC) was adopted to determine the deformation of the ground and the structure at the surface in three dimensions. Two high-resolution digital cameras whose relative position is very precisely known allow plotting of the observed 3D ground and structure surfaces for each of two speckle images taken simultaneously with the left and right camera of a stereo system. Then for two different speckle images, it is possible to determine the shape, the displacement, and the strain of the specimen's surface between different acquisition instants. In addition, this method provides an accurate result with an error smaller than 0.1 mm for a "discontinuous" surface such that of Fontainebleau sand.

The building model

A building model was created to investigate the effects of the soil-structure interactions during the formation of a subsidence trough. The geometry chosen for the building was inspired from the existing database of buildings damaged by mining subsidence in the east of France. A conventional house was considered as shown in Figure 8. This realistic but complex 3D prototype scale model was simplified for defining the small-scale model. In the first place, we considered only an elastic behavior and modeled an "equivalent slab". The procedure to simplify an individual house to a small-scale model is successively indicated on Figure 8. Three steps were carried out for the simplification procedure.

First, the equivalent slab was determined in such a way that the bending stiffness and the axial stiffness of this slab are closely equivalent to the 3D structure. We then intentionally chose to reduce the stiffness of the slab in both directions to exacerbate the strain in the structure. Both stiffness are approximately reduced by half. The last step is to transform this very simple structure into a small scale model by respecting the factors of the scaling laws. The structure model presented in Figure 8 is indeed a Usection slab made of polycarbonate, the interior part of which is composed of lead powder in plastic bags. This allows the model to present a stiffness and a stress transmitted to the ground approximately equivalent to those of the prototype. The 5 mm width of the edge is designed to be visible to the camera for measurement during the test. The main characteristics of the structure model and the equivalent prototype for a geometric scaling factor of 1/40 are shown in Table 2. The Table 2 also gives the value of relative bending and axial stiffness determined by equations 5 and 6 respectively.

The soil Young's modulus of 5 MPa was chosen because it corresponds to the low stress (< 1 kPa) applied on the soil below the structure. According to modification factors from the design curves of Potts & Addenbrooke (1997) as shown in Figure 5, the building model can be considered stiff for relative axial stiffness (M εhc and M εht close to zero) and quite flexible for relative bending stiffness in hogging zone (M DRhog ≈ 0.5) and very flexible in sagging zone (M DRhog > 0.5) for an eccentricity ratio e/B = 0.5 (position of structure center relevant to trough center).

To perform a test, the sand is placed manually in the physical model by layers of 0.15 m until the desired overburden height 0.3 m is reached. Each layer is compacted to obtain a certain level of density (a series of penetrometer tests was conducted to estimate the soil mass density, which is considered to be from average to dense). After a horizontal level of the ground surface is ensured, the small-scale building model is laid directly on the ground without any foundation system. The transfer of deformation to the building model is only done by the friction between ground-building interfaces, which limits the interaction between the soil and the structure. However, this concept makes it very simple to reproduce the tests. The image acquisition by the two cameras is then started with a frequency of 0.5 Hz. The motion program of the jack is launched at the same time a first snapshot is taken by the cameras. The jack is moved downwards at a rate of 0.125 mm/s for a total displacement of 30 mm. At the end of the test, the deformations of the ground and the structure are computed by Vic3D using the DIC method. It is important to note that the building model is removed from the ground surface to capture the displacement of the soil underneath the model at the final state of the test.

EXPERIMENTAL RESULTS

The aim is to discuss the soil-structure interaction and the building damage estimation using the existing methods (see section 3.2). Only one position of the structure relative to the subsidence trough will be analyzed herei.e. position 1, where the structure undergoes both compressive and tensile strains (see Figure 1). Position 1 was chosen because the structure was found to be the most affected by the ground movements (according to [START_REF] Caudron | A large 3D physical model: a tool to investigate the consequences of ground movements on the surface structure, 14[END_REF].

Greenfield ground movements

The ground movements in greenfield condition must be characterized and used as the reference for investigating the effect of soil-structure interaction.

Four identical tests were performed in order to ensure a good repeatability of the results. The variation of results from one test to another is not remarkable, but cannot be neglected. The configuration of the test is as shown in Figure 7, and the average global characteristics of the four tests are summarized in Table 3. As for the greenfield conditions, four tests were carried out for the configuration taking into account the building model. The result of one of the tests is presented in Figure 9, where the difference between the movements in greenfield and with building as well as the behaviour of the building can be seen. The results shown in Figure 9 are a slice from a 3D-surface of subsidence trough following a particular axis passing through the center of the building to that of the trough.

We can clearly see that the building causes important modifications in the soil displacement at ground level. The trough, formerly symmetrical under greenfield conditions, shows clearly some dissymmetry with the presence of the building.

Due to the load of the building and the contact with the ground, the trough size is enlarged. As a result, the ground slope is lessened. The building seems to reduce vertical ground displacement slightly. It is abnormal that the horizontal ground displacement in configuration with building is greater than that in greenfied. This could be explained by the fact that the building is located in the maximum slope area and that the ground could deform as the building model is removed to observe the subsidence trough.

The building behaviour deforms as a rigid body, since the horizontal displacement is almost a constant value along its length. Depending on the relative displacement of the ground in relation to the building model, three different zones may be distinguished.

In the center part of the trough, loss of support appears as the ground falls off the building due to a more substantial displacement. At the center part of the structure, the contact between the ground and the building is kept, and so important deformation of the building occurs, since the transfer of ground strains to the building is done by the friction from contact between the two surfaces. For the last part, the rigid body rotation of the building creates heave from the ground surface. The ground deformations in greenfield and with building, and the building deformation are compared and synthesized in Table 4.

Determining the deflection ratio of the structure is not straightforward. Two methods were then adopted to estimate the deflection ratio of the structure. First, it is determined by taking the maximum deflection (Δ) relative to the line connecting two reference points on either side (c,d) divided by the distance between the two reference points (see Figure 10). In another way, the deflection ratio of the structure is measured by the fraction of the maximum vertical distance (∆ h ) between the trough and the building model over the distance between two reference points (a,b) defining the contact zone. The range value obtained from both methods is presented in Table 4. The deflection of the ground is measured following the Potts & Addenbrooke approach (see Figure 4). The slope and strain present average values determined by differential vertical displacement and horizontal displacement respectively at the location of the building (c,d).

Figure 9. Displacements of the ground with and without structure at the surface and those of the structure (value at prototype scale) Figure 10. Building deformation within the hogging zone It can be seen that the measured deformations of the building are relatively small compared to those of the ground. Reasons could be that the building is just laid on the ground (no foundation system), the axial stiffness of the building is relatively important, and finally the building model is designed to perform in the elastic domain. The Potts & Addenbrooke (1997) approach shows that the ground strain can transfer from 0 to 10% to the surface building (see Figure 5b). If we consider the damage estimation chart of [START_REF] Wagner | Surface effects of total coal-seam extraction by Underground mining methods[END_REF], in which the ground strain (varies from 0 to 1%) is a damage estimation criteria, and Burland's chart (1995) based on structure strain (varies from 0 to 0.3%), we could estimate that, for a building 10 m long, the transfer from ground strain into the building is about 26%. Around 6.2% of the ground strain is transmitted to the building according to the results of physical modeling (Table 4)

The soil-structure interaction is thus clearly identified: the strains measured on the building are very different from those measured in greenfield condition, and that effect has to be taken into account for predicting building damage.

Building damage assessment

As shown in Table 5, with the data obtained from the physical model tests, the damage of the building model can be evaluated using the two-parameter charts presented in section 3.2.

The empirical chart of Wagner & Shümann (1991) (see Figure 2), which assumes that the building has no effect on the settlement prediction and the damage parameter (horizontal ground strain), is calculated using the predicted greenfield movements, makes it possible to estimate the risk level of the building model, which falls in severe to very severe damage. This is clearly an oversimplification of reality, as the stiffness of the building affects the transfer of ground strain to the surface building.

According to the analytical chart of Burland (1995) (see Figure 3) and with the movements of building models measured from the test, the likely damage of the building model is found from severe to very severe degree.

The Burland chart also makes it possible to predict the building damage from greenfield movements with the effect of building stiffness developed by Potts & Addenbrooke (1997). For each combination of the soilstructure relative bending and axial stiffness (see Table 2) and eccentricity ratio (e/B = 0.5), modification factors can be taken from the design curves of Figure 5. Then, the greenfield values of the deflection ratio and horizontal strain in hogging are multiplied by the respective modification factors (equations 2 and 4) to obtain those likely to be imposed on the structure. The combinations of hogging deflection ratio and tensile strain were then input into the Burland chart to quantify the likely building damage. This also estimates the damage of the structure from severe to very severe. Hence, the degree of building damage resulting from any estimation approach is from severe to very severe.

A good agreement between the building tensile strains determined from direct measurement and from greenfield movements with soil-structure interaction effect can be drawn from the Table 5. However, an important difference in building deflection ratios can be observed: about 4 times of the deflection ratio predicted from Greenfield incorporating soil-structure interaction according to Potts & Addenbrooke greater than that from the measurement. The cause could be explained by the hypothesis adopted in Potts & Addenbrooke approach that the building follows the ground movements, which is different from the behavior observed on the building model where the support loss occurs (see Figure 9). In addition, the elastic behavior of the building model can only be subjected to a maximum deflection ratio of 0.6%.

CONCLUSIONS AND PROSPECTS

The soil-structure interaction has been studied using 3D physical modeling making it possible to simulate the subsidence trough at the ground level. A simple structure model representing an individual house was designed to observe the importance of the soil-structure interaction to be taken into account for predicting structure deformation. The importance of building stiffness was proved by comparing the amount of strain transferred from the ground to the building obtained from experiments and from the design curves developed by Potts & Addenbrooke (1997).

The building damage has also been assessed using an empirical method applied for mining and another method used for tunneling. The latter makes it possible to estimate the damage based on the measured structure deformation and the modified greenfield movement with soil-structure interaction effect using the Potts & Addenbrooke approach. The degree of building damage was found to be severe to very severe for the different damage estimation methods. This shows a great coherence between them. However, it should be noted that the very critical damage obtained from the Burland chart is not due to the tensile strains, but to the excessive deflection ratio. The building strain determined from the greenfield taking the building stiffness effect into consideration is close to that from the experiments.

The physical model still needs to be developed. The present state of the model provides very significant deformation of the ground surface. Hence, several jacks should be installed to increase the variety of subsidence trough geometry and amplitude that the apparatus is able to achieve. In parallel, the building model should also be improved to more closely represent a real building with a foundation system and the upper part of the structure.

Table 1 .

 1 Relationship between category of damage and limiting tensile strain (after Boscardin & Cording 1989) Cording (1989) describe the damage corresponding

Figure 2 .

 2 Figure 2. Damage estimation according to[START_REF] Wagner | Surface effects of total coal-seam extraction by Underground mining methods[END_REF]

Figure 4 .

 4 Figure 4. Definition of deflection ratios (DR) in sagging and hogging (Potts & Addenbrooke 1997) Figure 5. Design curves for modification factors for (a) deflection ratio and (b) horizontal strain (Potts & Addenbrooke, 1997)

Figure 6 .

 6 Figure 6. Overall view of the 3D physical model

Figure 8 .

 8 Figure 8. Simplification procedure of a 3D individual house to an equivalent small-scale slab (dimensions not to scale)

Table 2 .

 2 Characteristics of building model and corresponding equivalent prototype

	Characteristics	Model	Scaling	Equivalent
			factor	prototype
	Width (m)	0.25	40	10
	Length (m)	0.25	40	10
	Total Height (m)	5.5E-3	40	0.22
	Young's modulus	2200	40	88000
	(MPa) Axial stiffness (MN)	0.66	40 3	42240
	Bending stiffness	2.03E-7 40 5	20.79
	(MN.m)
	Relative axial	1.06	-	1.06
	stiffness (-)
	Relative bending stiffness (m -1 )	1.67E-4 -	1.67E-4

Table 3 .

 3 Characteristics of greenfield ground movements (values at prototype scale)

	Characteristics	Notation Average
			value
	Max. subsidence (m)	S v,max	0.92
	Max. slope (%) Max. compressive strain (%) Max. tensile strain (%)	α max ε hc,max ε ht,max	23 -9 17
	Min. of curvature radius (m)	R ,min	6.5
	Distance from inflexion point to	i	4 to 5
	trough center (m)
	5.2 Soil-structure interaction

Table 4 .

 4 Comparison between the ground deformation determined in greenfield and with building, and the building deformation (values at prototype scale)

	Characteristics	Green-	With	Building
		field	building
	Max. Subsidence (m)	0.92	0.91	0.55
	Slope (%)	8.93	8.55	8.2
	Tensile Strain (%)	3.7	4.8	0.23
	Deflection ratio in	1.7	1.04	0.48 to
	hogging (%)			0.53

Table 5 .

 5 Damage estimation of building model

				Damage parameter
	Damage		L	ε ht	DR hog	Degree of
	estimation		(m)	(%)	(%)	damage
	method
	Wagner &		10	3.7	-	Severe to
	Shümann (1991)				very severe
	Burland (1995)		10	0.23 0.48 to	Severe to
	using measurements	1			0.53	very severe
	Burland (1995)	10	0.18 2.04	Severe to
	using Potts &				very severe
	Addenbrooke (1997) 2
	L: building length; ε ht : horizontal ground or structure strain;
	DR hog : deflection ratio of structure in hogging zone 1 structure deformation is obtained from measurements 2 The structure deformation based on the ground movements in
	greenfield and the soil-structure interaction effect