A Emonet

M A Bues

C Oltean

Bour

P Le Thiez

Migration of Trichloroethylene in a homogeneous porous medium. 2D numerical study

The infiltration of hydrocarbons and chlorinated solvents in soils and groundwater represents one of the most serious risks of long-term contamination. Transport of Dense Non-Aqueous Phase Liquids (DNAPL) and especially chlorinated solvents is controlled by their physical and chemical properties. The aim of this work is to better understand the behaviour of Trichloroethylene (TCE) in a homogeneous porous medium by using both experimental data and numerical simulations. The study focuses on the influence of capillary fringe and differential water saturation on the behaviour of TCE phases. The experimental data of capillary pressure curve versus water saturation obtained for a homogeneous porous medium is used in the numerical study carried out with SIMUSCOPP package and is compared with simulations performed without taking into account capillary fringe.

INTRODUCTION

Dense Non Aqueous Phase Liquids (DNAPL) encountered in subsurface can have a variable physical and chemical behaviour depending on the type of fluid present in soils (oil, water and gas) as well as the soil matrix [1]. It can be found under different chemical species: dissolved compound, organic liquid, volatile compound or adsorbed compound. Each form represents a pollution risk for groundwater and soil. A better understanding of these mechanisms will allow us to predict the change of front propagation and help us to elaborate remediation processes.

Mass transfers between DNAPL and water, during the dissolution process, have been studied by authors during the 1990's ( [START_REF] Miller | Dissolution of trapped nonaqueous phase liquids : mass transfer characteristics[END_REF], [START_REF] Powers | Theoretical study of the significance of non-equilibrium dissolution of non-aqueous phase liquids in subsurface systems[END_REF], [START_REF] Powers | An experimental investigation of non-aqueous phase liquid dissolution in saturated subsurface systems : steady state mass transfer rates[END_REF], [START_REF] Geller | Mass transfer from non-aqueous phase organic liquids in water-saturated porous media[END_REF], [START_REF] Powers | An experimental investigation of non-aqueous phase liquid dissolution in saturated subsurface systems : transient mass transfer rates[END_REF]). Dissolution processes take place in several stages. The first is the movement of a compound inside the NAPL phase towards the fluid interface. In this case, chemical transfer occurs at the interface between fluids. Diffusive processes will cause mobilisation of dissolved compounds, away from the interface.

The maximum concentration of dissolved compound might be lower or equal to solubility but this process may be limited by different causes ( [START_REF] Abriola | Modeling multiphase migration of organic chemicals in groundwater systems. A review and assessment[END_REF], [START_REF] Powers | Theoretical study of the significance of non-equilibrium dissolution of non-aqueous phase liquids in subsurface systems[END_REF], [START_REF] Powers | An experimental investigation of non-aqueous phase liquid dissolution in saturated subsurface systems : steady state mass transfer rates[END_REF]) such as: non uniform flow of aqueous and organic phases due to aquifer heterogeneities, influence of equilibrium concentration on the composition of aNAPL mixture, etc.

Mass transfers during volatilisation have also been studied ( [START_REF] Sleep | Modelling the transport of volatile organics in variably saturated media[END_REF], [START_REF] Griffol | Contaminant migration in the unsaturated soil zone: the effect of rainfall and evapotranspiration[END_REF]). Sleep and Sykes [START_REF] Sleep | Modelling the transport of volatile organics in variably saturated media[END_REF] showed that movement of volatile compound due to advection is influenced by density driven flow. Cho et al. [START_REF] Cho | Simulating the volatilisation of solvents in unsaturated soils during laboratory and field infiltration experiments[END_REF] and Jellali et al. [11] have investigated the dissolution process of TCE vapour after water infiltration.

In our study, we focus on the influence of water saturation on mass transfers between TCE organic phase and dissolved or volatile phases during infiltration across the unsaturated zone up to reaching the saturated zone. Li this study, we do not take adsorption into account. The behaviour of TCE depends on water saturation i.e.: unsaturated zone, capillary fringe or saturated zone. Thus, it is important to precisely evaluate the local water saturation of the porous medium. The determination of capillary pressure data provides us with relative permeabilities which control TCE behaviour in the capillary fringe.

Several theoretical models based on capillary pressure, water saturation and relative permeabilities will be introduced first in this paper. Subsequently, it is explained how the experimental retention curve (capillary pressure versus water saturation) was obtained in order to determine the water saturation change along the capillary fringe. In the final section, we compare SIMUSCOPP simulations implemented with and without considering capillary fringe.

GOVERNING EQUATIONS FOR THE OBTAINING OF CAPILLARY PRESSURE AND RELATIVE PERMEABILITIES DATA

Capillary pressure and relative permeabilities are important data in multiphase fluid transport in soils but are quite difficult to obtain experimentally. Therefore, several theoretical models for predicting capillary pressure and relative permeabilities have been proposed. In this study, we retain the models proposed by Van Genuchten -Mualem and Parker.

The experimental soil-water retention curve of Van Genuchten (1980) [START_REF] Van Genuchten | A closed-form equation for predicting the hydraulic conductivity of unsaturated soils[END_REF] yields the parameters necessary to predict relative permeabilities using Mualem's model. The general equation (1) links effective water saturation S* (Eq. 2) to pressure head (h, h>0): where n, m and a are curve shape parameters determined with the soil-water retention curve.

Van Genuchten showed that P is the inflection point and for P we have also:

s.- {JSsaturated -Oresidualj (5) dQogh) ifO<Sp<\ thenm = l-exp(-0.8S£) and if Sp>\ then m = 1-- +4 + ^. Sp Sp 2 Sp 3
Then relative permeabilities for the couple water-air for a soil are given by equation ( 6):

) = S w * uz [l-(l-S w * Um ) m Y with0<m<l andEq. 3 (6)
Equation ( 6) is valid for two-phase, air-water system. This can be adapted for a threephase system i.e. air, water and NAPL [START_REF] Parker | A parametric model for constitutive properties governing multiphase flow in porous media[END_REF]. In such a system, NAPL is supposed to have an intermediate wettability between air and water. We assume in this paper that this formulation is applicable for all NAPL. For each couple the two-phase formulation for the different fluids in a three-phase system yields:

S *f = [1 + (Ughf,)"T m if hy > 0 and S • * = \iîh ij <0
for i, j = a, o, w with i £j (7)

S *r = [1 + (aaokofr if K, > 0 and S *T =lifh ao <0 (8)
where the total liquid saturation in a three phase liquid system is defined by S*t and S *'j is expressed as:

J 1 C 1 O 1 -Oresidual I -Oresidual
Other studies have shown that NAPL relative permeabilities in a three-phase system depend mainly on gas and water saturation. According to Parker [START_REF] Parker | A parametric model for constitutive properties governing multiphase flow in porous media[END_REF], if the air-water relative permeabilities are given by equation ( 6), then, for the other fluids, these equations can be expressed as: 2m with C a gas slippage correction coefficient (11)

ka = CSa * 1/2 (1 -S, * llm )

{lkro = (S* t -S*J n {l-S^1"]" -[l-S*]'" 1 ] 1 } (12)

These models enable us to determine relative permeabilities for the three phases from capillary pressure and water saturation data. The data obtained experimentally and processed by this method is implemented in SDVIUSCOPP. Its influence will be shown in paragraph 4.

DETERMINING THE RETENTION CURVE

The porous medium used in the whole experimental study is a natural sand from the Moselle river, eastern France. It is mainly composed of quartz with a small proportion of feldspar. Its main characteristics are as follows: effective porosity of 0.38; hydraulic conductivity close to 2.1(rtn/s and mean grain-size between 0.4 and 0.6 mm.

The experimental apparatus (Figure 1) used for the determination of the retention curve consisted of a stack of 8 plexiglas cylinders (9 cm in diameter and 7 cm in height). The assembled column height is about 55 cm. It is first filled with water then sand is progressively added in 400 g portions. The column is struck several times with a hammer in order to optimise grain sand distribution and to avoid compressive effects. Tensiometer probes are attached on the cylinder; one end is buried the porous medium and the other end is linked to a pressure gauge. The change in pressure inside the porous medium during the experiment is recorded. Three tests have been conducted in the same experimental conditions in order to test the retention curve. The results of these tests are presented in Figure 2. The variation of each parameter is between 15 % to 18 % for Smin and n and around 6.5% for m and a. Values obtained for test 2 are intermediate compared with test 1 and test 3, so we may consider that it is representative of the other tests.

These parameters enable us to estimate the relative permeabilities for each two-phase fluidcouple in a three-phase system. In SDVIUSCOPP, hysteresis is not considered. Therefore, we will focus our attention on the drainage curve and not on that of imbibition. These experimental and analytical values have been implemented in SDVIUSCOPP and compared to a simulation which does not consider capillary pressure.

SIMULATION OF TCE MIGRATION WITH SIMUSCOPP

SIMUSCOPP is a numerical modelling software package developed by the French Institute of Petroleum (D?P). It solves a set of coupled transport equations for fluid velocities, pressure and concentration of chemical species in multiphase flow. Several chemical species can be found in one, two or three fluid phases. The governing equations are based on mass balance and momentum equations, except that those of Darcy are used in place of those of Navier-Stokes. Two cases can be considered, whether or not mass transfer between phases and fluids is taken into account. The geometry considered in these simulations uses regular radial coordinates.

The set-up simulated represents a column, 20 cm in internal diameter and 100 cm in length. It is divided into 400 elements of 0.5 cm in radius and 1 cm in depth. We consider the column as a cylinder. Numerical simulations are applied to the angular width of each cell. Graphic representations are illustrated using a vertical section of the column where we Unsaturaicd Zone _o,3-Capillary Fringe consider only a radius. In this way, the injection of TCE which takes place in the centre of the column is represented in one side of the column. The characteristics of the porous medium are those indicated in paragraph 3. In this case, SIMUSCOPP is used in a multiphase structure; therefore relative permeabilities and capillary pressure data for each couple of fluids are required.

For each simulation, the unsaturated zone is situated at depths between 0 and 0.5 m, the capillary fringe between 0.5 and 0.8 m and the saturated zone between 0.8 and 1 m. 100 mL of TCE are injected over 2 hours in the upper part of the model. The injection zone (1.5 * 1 -cm 2 , i.e. 3 cells in direction X and 1 cell in direction Z) is situated 5 cm below the ground surface, in order to avoid volatilisation of TCE. The different species of TCE (organic phase, volatile and dissolved) are illustrated in Figure 4.

Five simulations experiments have been conducted: one without the implementation of capillary pressure data (test 0) and the others with capillary pressure and relative permeabilities obtained experimentally (test 1 to 3), as described in paragraph 2 and 3. The influence of each experimental capillary pressure -water saturation curve on water saturation distribution along the column is illustrated in Figure 3. In a simulation where capillary pressure is ignored (test 0), a distinct separation surface is obtained between the two zones: saturated and unsaturated, without any gradation in water saturation. The capillary fringe is absent. On the other hand, the use of the experimental curves induces the appearance of a capillary fringe whose height varies between 20 (test 3) and 30 cm (test 1). Moreover, as numerical simulations give similar results for test 1 and test 2, we can consider one of these tests (e.g. test 2) as the most representative. The influence of the experimental capillary pressure -water saturation curve on the behaviour of different species of TCE will be discussed for test 2. Numerical results of this simulation are illustrated and compared with those of test 0, in Figure 4.

Unsaturateâ Zone Capillary Fringe o.i- -02- -0.3- -0.4- -0.3- • 0.6 - --o.i- --0.2- --0.3- --0.4- --0.5- --0.6- --0.1- --0.2- --0.3- --0.4- --0.5- --0.6- --o.i- --0.2- --0.3- --0.
In Figure 4a, one can observe the evolution of oil saturation in time for test 0 (dark) and test 2 (light). For test 0, TCE penetrates rapidly into the porous medium and spreads through the lower part of the column (saturated zone). The migration is primarily vertical. Conversely, TCE in test 2 crosses the unsaturated zone more slowly with a more extensive lateral extension. When TCE encounters the capillary fringe it slows down and accumulates in the upper part of the saturated zone.

In Figure 4b and Figure 4c, the migration of volatile TCE and dissolved TCE is illustrated. These chemical species are represented according to its maximum molar fraction for TCE vapour (equal to 0.077) and maximum concentration in water for dissolved TCE (equal to TCE solubility, Le. 1100 mg/L).

Concerning the migration of TCE vapours, numerical simulations show markedly different behaviour between test 0 and test 2. For the first one, we observe that during the first hour, TCE vapours are located at the place of the TCE path, towards the centre of the column. Vapours then begin to spread laterally and after 3 hours they occupy almost the whole of the unsaturated zone. For test 2, in the same time interval, TCE vapours follow the slow migration of TCE in the column. Lateral extension is more extensive than in test 0 and vertical migration is slower. After 2 hours, the unsaturated zone is almost completely saturated with TCE vapours.

Figure 4c shows the evolution of dissolved TCE. For test 0, dissolved TCE is only located in the saturated zone during the first hour. After 2 hours of simulation, dissolved TCE is also found in the unsaturated zone. This is probably due to solubilization of vapours in contact with irreducible water in the unsaturated zone. In the saturated zone, dissolved TCE is present at the bottom of the column. For test 2, dissolved TCE, as for TCE vapours, follows the migration of TCE organic phase. The lateral migration is significant and the vertical migration stops at the upper part of the saturated zone.

CONCLUSION

In this study, we have shown the experimental results of determining capillary pressurewater saturation curve for a porous medium. For a homogeneous porous medium, similar experimental curves of capillary pressure versus water saturation have been obtained. In a numerical model, this experimental data and estimated parameters enable us to numerically study the migration of TCE with or without taking the capillary fringe into account. The behaviour of TCE in these two cases is markedly different. The presence of a capillary fringe causes a slowdown of TCE migration but a more extensive lateral extension and an increase in the velocity of contamination by TCE vapours and dissolved TCE in the unsaturated zone.

These numerical results will be compared with new experimental results concerning TCE migration in a 2D column.

Figure 1 : 2 :

 12 Figure 1: Experimental set-up.

Figure 3 :

 3 Figure 3: Distribution of water saturation in the column after simulation with SIMUSCOPP.

Figure 4 :

 4 Evolution of TCE migration versus time, for its different phases.

Table 1 :

 1

	Table 1 Test	S max	S min	m	n	a [m 1 ]
	Testl	1	0.1037	0.658	2.92	0.1032
	Test 2	1	0.0905	0.684	3.16	0.0919
	Test 3	1	0.0712	0.745	3.92	0.0930
	Mean	1	0.0885	0.696	3.33	0.0960


Abbass Fa

Y Ad

Christophe Didier
email: christophe.didier@ineris.fr

Modelling propagation of sinkhole, in both slow and dynamic modes, using the UDEC computer code. (Methodology and application)

Keywords: Sinkhole, structure, propagation, influence, modelling, static calculation, dynamic calculation, forecasting methodology

The propagation of sinkhole phenomena towards the surface is a fairly complex phenomenon, with the potential to cause damage to surface structures within the affected area. Using numerical modelling a methodology has been developed which makes it possible more closely to specify the conditions in which sinkhole forms and to propose a prediction model. The UDEC code is used. An actual case of sinkhole in the Paris region (Vachat, 1982) has been used as a basis for the study. We simulated the predicted collapse mechanism and the profile of the sinkhole itself. The influence of a surface structure was also examined. The modelling process in static and dynamic modes was developed in order to simulate the ways in which sinkhole forms. As regards dynamic modelling, we determine the relevant parameters: critical damping and fundamental frequency. The dominant point in this method is that induced fractures in the roof strata overlying underground excavations -the root cause of sinkhole -are created by hand. The fractures were induced using the criterion whereby a zone is formed in which the failure criterion is exceeded. The comparison between observation and the modelling results is very satisfactory. It may be noted that the numerical modelling and the methodology used enabled us to determine the shape and parameters of the roof "bell" and of the zone of sinkhole. We were also able to confirm the influence of external loads that might facilitate development of the sinkhole. A feature of dynamic (rapid) appearance of sinkhole is apparently that it takes place on a bigger scale. Incidentally such an approach could be used to examine the consequences for structures.

INTRODUCTION

When ground -frequently stratified -collapses owing to the presence of underground voids, particularly in shallow mines and quarries, the phenomenon propagates upwards in a bell-shaped profile and can induce sinkhole at the surface (fig. 1). Underground the bell-shaped profile is usually characterised by its height (h), its diameter (D) and an angle that serves to indicate the shape of the bell (a). At the surface, sinkhole is characterised by its depth (P) and diameter (<j)).

Observations of different cases of sinkhole have been used to propose empirical prediction equations (Vachat, 1982). Although pragmatic, this approach is not accurate enough and is limited to areas that are geologically uniform and hence fairly limited in extent. Approaches based upon the strength of materials are often used for these predictions [START_REF] Didier | Comportement des ouvrages historiques soumis à des sollicitations sismiques[END_REF].

These types of sinkhole are frequently the cause of more or less substantial damage depending on their dimensions and on the type of structure (fig. 1). Any attempt to improve the behaviour of structures, once sinkhole has appeared, is a tricky business. This article seeks to improve the prediction of sinkhole development and its consequences for structures (fig. 1) using modern methods (numerical modelling using distinct elements, UDEC). This new approach makes it possible to determine the probable profile of die sinkhole and to specify the influence of structures according to the behaviour of the ground.

Bell shape

Consequences 2 Qusugy

CHOICE OF METHODOLOGY PROGRAM AND

The failure mechanism involved in sinkhole, in the case of stratified rock, is die yielding of the different roof strata [START_REF] Peng | Predicting surface sinkhole for damage prévention[END_REF]. It is possible, using the finite element method in the UDEC code (Itasca), to localise the failure zones where the criterion is exceeded.

On the other hand it is not possible with UDEC to create induced fractures within diese zones. For simulating the failure of the strata, the possible fall of blocks, and the appearance of sinkhole at the surface, we propose to generate these fractures by hand.

We assumed that rock failure was likely to occur when plasticity developed along the cross-section of a given block. On the other hand if the areas of plasticity are discontinuous, we shall speak of partial failure (partial fracturing), and the block concerned is regarded as not undergoing any geometrical change. This is a progressive approach: we begin the process with the first block in the roof and then apply it to the others. The analysis terminates when no more new fractures can be created using the above criterion.

DESCRIPTION INVESTIGATED OF THE CASE

Description of the site and properties

For the purposes of this methodological study, we turned to an actual case of sinkhole in the Paris region. The quarry involved was shallow (8 m) and had been excavated using the room and pillar method. The depth of the block worked did not exceed 2 m, giving a ratio of overburden to gallery height of less than four. The sinkhole had a diameter of 2.5 m and a depth of 2 m (Vachat, 1982). The bell profile formed was characterised by an angle of 80° at the base. The overburden consisted of a series of blocks of limestone, marl and clay.

Simulating the case using UDEC

The aim is to produce a model to simulate the development of sinkhole. Sinkhole usually has a conical profile. The model adopted in two dimensions is based upon the assumption of plane deformations. This assumption is more pessimistic than that of axial symmetry (not available in UDEC), a fact confirmed by Mandel (1966) who compared the two solutions for a cylinder and a sphere. The equations for the stresses (OQ and CT r ) given below are for a semi-infinite cylinder of radius (a):

O t tunnel = G 0 (l-a 2 /r 2 ), 2?/1 7 ). tunnel = (1) (2)
The equations for the distribution of the stresses (Oe and a r ) around a sphere of radius (a) that represents the model with symmetry of revolution are:

r sphere = C o (l-a V r 3 ),
Figure 2 shows a comparison of the solutions for the two approaches. It can be seen that the difference is significant, particularly at a distance of one to two times the radius of the cavity Figure 2 : Distribution of stresses, for a = lm and Ob = -IMpa.

Law governing the behaviour of materials

The model is elastic-plastic with behaviour of the softening type (negative work-hardening). This represents the behaviour of the rock under stresses that are induced in the vicinity of the underground workings (Piguet 1983 and[START_REF] Fine | Le soutènement des galeries minières[END_REF]. In the plastic zone, the work-hardening variable depends on the plastic distortion (Prévost and Hoeg, 1975). The variation in cohesion (C), the angle of friction ((p), the tensile strength (Gi), and expansion (l|/) as a function of the plastic deformation (£ p ), are shown on figure 3. In certain cases, the initial values of C and (p fall by half with the value G r falling to zero, for a value of plastic deformation (8 P ) of 1% [START_REF] Tharp | Design against collapse of karst caverns[END_REF]. These factors are taken into account in our study.

Law governing the behaviour of junctions

Force-displacement equations are used in UDEC to indicate the behaviour of junctions. We used an elasticplastic law with variable stiffness (a function of the stress normal to the junction plane G n , jK = aG n b

). The values of (a) and (b) are functions of the stiffness (normal or tangential) respectively. They are given by the following equations:

=100<T jK a =100<T 2 n , ( 6
) (4)
Figure 4 shows the model adopted. It is recalled that the angle of friction of the junctions (j(f>) were obtained by Filliat (1984) and [START_REF] Fine | Le soutènement des galeries minières[END_REF]. We also assumed that the simple tensile strength (j<J t ) was zero, cohesion (jC) negligible, and the angle of expansion of the junctions as given in the UDEC manual [START_REF] Itasca | UDEC (universal Distinct Element Code). Version 3, Volume I et II. Mandel[END_REF].

For the purposes of this modelling process, it is necessary to estimate the tensor of the initial stresses (no measured data are available). We assumed that the ground was subjected to its own weight and therefore to a vertical stress of pyh and a horizontal stress of 0.3 pyh.

Results of the simulation (formation of sinkhole)

Figure 5 represents the potential zones of failure following simulation of the excavation of the gallery. It can be seen that the first block is liable to fracture in two zones. These failure zones are replaced by two fractures forming a central block and two parts that appear to overhang (figure 6).

10m

Need fo cf-& They conic t iï i ^ffe • y* t • i- ", 0 roof; -o < v B.ss i ori.
--Nwlrc/I oe/s J j Quarry I Figure 6 shows the configuration obtained, which will be modelled in its turn. This simulation is one of the intermediate stages moving towards final equilibrium.

After creating these fractures, a new stage of calculation is undertaken with a further analysis of the zones of plasticity in this intermediate stage. The presence of new zones means that the creation of new fractures is repeated up to equilibrium and the complete stability of the roof with or without the appearance of sinkhole at the surface (fig. 7).

Figure 8 shows the final shape of the cavity after the seven phases of calculations. Equilibrium has been reached following the failure of the 7 th block (fig. 8). It will be noted that the failure zones cannot generate further blocks.

In conclusion, no sinkhole has reached the surface. HoeScNo.9

Block No,a 10m 7,6m In this configuration, the bell shape is characterised by a height (h) of 6 m, a base angle (OC) of 76° and a diameter at the base (D) oi l.6 m.

Table 3 shows the lengths of different blocks in the roof.

It can be seen that the span of the blocks diminishes as we approach the surface.

Block

Simulation of a surface structure in the model

To examine the influence of aggravating factors on the speed of propagation of die sinkhole to the surface, and on the resulting consequences for the structures affected, we simulated the presence of a four-storey building above the underground cavity. Figure 9 shows a plan view of one floor of the building with the data necessary for calculating the loads acting on the pillars Pi, P2 and P3. The transverse section (X-X) shows the loads transmitted to the ground at the foundations (A, B and C). The transmitted loads are calculated in a simplified manner. The results show a load of 100 kN under the foundations A and C and a load of 200 kN under the foundation B.

Analysis of the final configuration

The procedures indicated earlier are applied to the new model, to establish the same criterion for generating fractures. Figure 11 shows the configuration of the sinkhole. This is generated beneath foundation A in the centre of the model.

The bell-shaped sinkhole formed is characterised by: a height (h) of 8 m, a base angle (a) of 76° and a diameter at the base (D) of 7.6 m. The form of the sinkhole is characterised by a diameter at the surface (<))) under foundation (A) of 1.5 m. The mean lengths of fractured spans in the roof blocks (with structure) are shown in table 4, from the bottom up. Comparison of tables 3 and 4 shows that the blocks of depth gi'eater than 3m are unaffected by the presence of the building. The lengths of the failed spans are the same in the two models. This may be explained by the width of the foundations (1=1 m), and the influence of the loads does not exceed a depth of 3 times 1 (= 3 m in the particular case examined). For the other blocks however, the influence of the building shows up clearly as an increase in the length of the failed spans and by the occurrence of sinkhole (fig. 11).

SIMULATION IN DYNAMIC MODE

The aim of this simulation is to investigate the rapid formation of sinkhole. We assumed that there is not sufficient time for the rock mass to damp out the shocks induced by the failure of the roof blocks. This type of calculation requires data such as the fundamental frequency (f p ) of the model and its critical damping (^c).

Sizing the dynamic model and the calculation mesh

In order to meet the requirements of the dynamic mode calculation, the model has to be resized: the dynamic calculation is affected to a greater extent by the size of the model and that of the mesh (due to wave reflection).

Figure 12 shows the cross-section of the model together with the proposed mesh (table 5) which will be validated numerically after calculating f p .

For the dynamic mode calculation to be reliable, the maximum size (T) of the zones (or elements) should not exceed approximately 1/8 of the wavelength (À.) to be

Calculation of critical damping (£ c )

The value of critical damping can be selected from experimental values found in the literature. The value taken was between 2.5 and 5%. For f p < 10 Hz, it is necessary to multiply t,c by two (El-Shabrawy 1996).

Calculation of the fundamental frequency (f p )

Calculating the value of f p involves looking for the fundamental frequency of the system (the lowest frequency). To do this, we carried out a calculation without damping by applying only static loading, and then recorded the system response (speed or displacement) at any point (fig. 13). Since there is no damping, the system oscillates freely and its fundamental frequency can be determined from observation of its oscillations (El-Shabrawy, 1996).

Zero damping should always be introduced after simulating the excavation, however one may wonder whether this phenomenon takes place before or after static equilibrium. The dynamic effects of an earthquake are often introduced after static equilibrium. In our case, zero damping is applied with respect to two options: the excavation carried out without induced fractures (fig. 5) and the second after formation of the bell shape (fig. 8).

Calculation of f p was done at two points in the model P (10; 9.5) and Pi(0; -45) as shown in figure 12. The point P is placed above the gallery in order to optimise f p in the zone of high perturbation. Pi however is placed at the bottom of the model in order to optimise f p in the zone of low perturbation. Figure 13 shows how the velocity varies as a function of time at the points P and Pi for zero damping. The two figures (13-a) and (13-b) are obtained for the first option. The resulting values for f p are 17 and 20 Hz respectively. The two figures (13-c and 13-d) correspond to the second option. The result is values of f p of 4 and 10 Hz respectively. We can see therefore that f p varies between 4 and 20 Hz (fig. 13). This variation is due to the method of applying zero damping and to the behaviour of each stratum (block), particularly in the zones that are affected or perturbed by the excavation.

, I! f

M IN il sis in mm ! pi! I f E =17Hz 13-a 13-b f p =10Hz
13-c 13-d Figure 13: Calculation of f p using the method of Ç c -0.

The mesh used is in fact verified for the minimum velocity. From figure 13 it can be seen that the minimum velocity of the wave is that of figure 13-d). The time the wave takes to reach the point Pi is 0.1s. The distance between the origin of the frame of reference (fig. 12) and Pi is easily calculated (d = 45 m). We conclude that the velocity of the wave is about 450 m/s. We can say then, from equation ( 10) and putting f max ( = fP = 10 Hz) that the mesh size should be < 5.6 m (T < 1/8.450/10 => T < 5.6 m). Accordingly the mesh used (table 5) meets the conditions required for a dynamic calculation.

Calculation method

The dynamic simulation was done for four modes, each mode being characterised by the values of < f y and f p (fpnmx-20 and fpmm = 4). It will be noted that modes 1 and 2 correspond to dynamic behaviour only after the appearance of the sinkhole, while modes 3 and 4 correspond to dynamic simulation of sinkhole and the excavation. Table 6 shows the parameters used for each mode.

Simulation

Analysis of the final configuration

We applied the procedures indicated in paragraph 5. The results for the different modes of dynamic simulation (table 6) concern the shape of the sinkhole "bell" and are given in table Table 7 and figure 14 show the shape of the sinkhole bells together with the relevant parameters. Depending on how the model is simulated dynamically, the diameter of sinkhole at ground level varies between 1.4 and 3.2 m.

The realistic cases (modes 3 and 4) correspond to bigger sinkhole bells than cases (modes 1 and 2). It may be noted that the results (sinkhole bells) as well as the surface movements vary substantially according to the simulation method used (pessimistic or optimistic) and the critical damping value applied.

Analysis of the final configuration with structure

The configuration of a void with a structure is modelled in dynamic mode ( § 7.5). Table 8 and figure 15 show the results: it can be seen that the sinkhole at the surface has a diameter varying from 2.2 to 3.2 m. The influence of the structure appears in modes 1 and 2 corresponding f p =4 and low critical damping. The simulation of the model for dynamic excavation (modes 3 and 4), despite the increase in £, caused considerable perturbation in the roof blocks and very high stresses. The influence of the loads in this type of structure on the block failure processes in negligible. The influence of the structure on the failed spans is also negligible.

COMPARISON OF THE SIMULATION RESULTS WITH OBSERVATIONS

We shall refer to the case investigated by Malakoff: this was sinkhole 2.5 m in diameter at the surface and 0.5 m in depth. In order to compare the results obtained with the model and observations, we show table 9 which gives the values obtained for the shape of the bell and sinkhole for all the calculation modes (mean value in dynamic behaviour):

Comparison with the calculation in static mode

For static behaviour, and contrary to observation, no sinkhole was observed when no structure was present (table 9). The presence of a structural load led to the formation of sinkhole at the surface, but of smaller diameter (1.4 m rather than 2.5 m). The gap between observation and simulation is narrower. This confirms the influence of external loads on the formation of sinkhole. The values of 0C (=76°) and D (=7.6 m) can be used to calculate the maximum geometrical height of the "bell" (li c ) [given by Vachat (1982)]. We have h c = D.TgOC/4, or 7.6 m. This value is the order of magnitude of those obtained by simulation. The following remarks are relevant: the ingress (or presence) of water, the loads applied and collapse events at the surface all play a predominant role in the ascent of the "bell" (Vachat, 1982). However, these factors were not modelled.

Comparison with the dynamic calculation

The results in dynamic mode are closer to observation. The dimensions of the sinkhole are greater. These dimensions are increased further by the combined effects of the structure and the speed of appearance of the sinkhole. The diameter is higher than or equal to that observed in site.

INDUCED CONSEQUENCES (OR GROUND-STRUCTURE INTERACTION)

In order to evaluate surface movements and the consequences due to the appearance of sinkhole, it is necessary to introduce the parameters of the subsidence depression (displacement, deformation, gradient, and so on), pointing out that the appearance of sinkhole at the surface causes a change to this depression and, naturally, to these parameters [START_REF] Bibliography Al-Heib | Effets et conséquences des mouvements du sous-sol sur le bâti. Synthèse des principales méthodes de prévention et de protection[END_REF]. These changes may be limited to the vicinity of the sinkhole diameter (discontinuous movements). Beyond this, the characteristics of the depression and the rate of change of its parameters normally remain unchanged (fig. 16).

Figure 16 shows the half-profiles of vertical displacement (d 7 ) and horizontal displacement (d x ). It will be seen that the slopes of the vertical and horizontal displacement curves in dynamic mode are greater than those in static mode. Table 10 shows the displacements and deformations at the locations of the foundations A, B and C: It can be seen that:

^ the dynamic calculation causes heave at foundation (C), apparently due to substantial collapse of the ground after appearance of the sinkhole;

^ the maximum values induced at foundation C are due to the static calculation;

S the deformation induced at the foundations is practically negligible. This is confirmed in practice (rigid foundations prevent or minimise the deformation of the ground, see Deck 2002). These remarks suggest that the dynamic behaviour induces more substantial consequences.

In table 10 the term sinkhole indicates that sinkhole has completely damaged foundation (A). Its stability could be assured by appropriate methods [START_REF] Bibliography Al-Heib | Effets et conséquences des mouvements du sous-sol sur le bâti. Synthèse des principales méthodes de prévention et de protection[END_REF], for example:

S filling the underground cavity by injection;

S building a special foundation to prevent this occurrence (fig. 17a);

• S building the structure over an overall raft (fig. 17b) or with a gantry structure (fig. 17-c) to resist the collapse of the ground under the foundation. As far as foundations B and C are concerned, their size will accommodate the relevant movements.

CONCLUSIONS AND OUTLOOK

The study was able to:

S simulate and provide understanding of the phenomenon of the ascent of the "bell" as a result of the manual creation of induced fractures;

S predict the form of the propagation of the sinkhole towards the surface;

S show that the presence of a structure at the surface can accelerate the movement of the void towards the surface (aggravating factor);

S confirm that the sudden failure of roof blocks applied in the dynamic calculation generates resonance forces that enhance the growth of the underground void;

S predict the shape of the settlement depression and deduce the movements induced in the affected structures.

The calculation procedures used represent a novel method for predicting the appearance of sinkhole and for evaluating its induced consequences in affected structures. This modelling approach could not be generalised unless applied to several cases. The aim is dierefore to extend this methodology so as also to encompass parameters related to the environment (ground, structures, and so on).

Figure 1 :

 1 Figure 1 : Development of sinkhole and induced consequences.

Figure 3 :

 3 Figure 3 : Relations between C, (p, G t , \|/ and S p .

Figure 5 :

 5 Figure 5: Location and creation of induced fractures.

Figure 6 :

 6 Figure 6: Intermediate phase of creating induced fractures.

Figure 7 :

 7 Figure 7: Intermediate phase of the final calculation.

Figure 8 :

 8 Figure 8: Final configuration for the case examined.

Figure 9 :

 9 Figure 9: Illustration of the simulated building.Owing to the difficulty of simulating structures in UDEC, the load of the structure is simulated by applying normal stresses, like boundary conditions, equivalent to the loads calculated in the previous paragraph. These are placed (fig.10) respectively above the centre of the gallery (100 kN), above the right-hand boundary (200 kN) and 4 m from the right-hand boundary (100 kN), each 1 metre wide.

Figure 10 :

 10 Figure 10: Incorporating the building in UDEC.

Figure 12 :

 12 Figure12: Cross-section and mesh used for dynamic calculation.

Figure 14 :

 14 Figure 14: Final dynamic configuration.

Figure 15 :

 15 Figure 15: Final configuration in dynamic mode with structure.

Figure 16 :

 16 Figure 16: Variations in the half-profiles of displacements.

Foundations

  Parameters for deformations induced by the sinkhole.

FigureFigure 17 :

 17 Figure 17-a: Special foundations

Table 1

 1

	shows the

Table 1 :

 1 Mechanical and geometrical properties of materials.

	Junction No	1	2	3	4
	maxjk n (Mpa)	30	40	60	80
	maxjk s (Mpa)	15	20	30	40
	minjk n (Mpa)	20	30	50	60
	minjk s (Mpa)	10	10	30	30
	J<P (°)	20	20	20	22
	jo, (MPa)	1	1	1	1
	jC (MPa)	100	100	100	100
	V(°)	10	10	10	10

' ' • i . ' . " ÎX 1 -é r ' W«..r .. UlHi~~ . . sfy.% 7%*çT-Figure 4: Cross-section and materials in UDEC. The properties of the junctions between strata are shown in table 2. The maximum and minimum values of normal (maxjKn and minjK n ) and tangential (maxjK s and minjK s ) stiffnesses are functions of the values of K and G for the corresponding materials (based upon experience).

Table 2 :

 2 Mechanical properties of junctions.

Table 3 :

 3 Lengths of failed spans.

	No.	Block depth (m)	Length of failed span (m)
	1	7.25	7.60
	2	6.00	7.10
	3	5.00	6.20
	4	4.20	5.40
	5	3.00	3.84
	6	2.50	2.50
	1	2.00	1.00
	8	1.00	0.00
	9	0.00	0.00

Table 4 :

 4 Lengths of failed spans (with structure).

			1,6m	c
	10m
		_ {	Quarry	|
	Figure 11: Final configuration with structure.
	Block No.	Block depth (m)	Length of failed span (m)
	1		7.25	7.60
	2		6.00	7.10
	3		5.00	6.20
	4		4.20	5.40
	5		3.00	4.84
	6		2.50	4.30
	7		2.00	4.10
	8		1.00	2.80
	9		0.00	1.50

Table 7 :

 7 7: Shape of the bell and the sinkhole in dynamic mode.

			Simulation	Sinkhole 4>(m)	h(m)	Shape of bell <x(°)	D(m)
				Mod.l	1.4	8		76	7.6
				Mod.2	1.0	8		76	7.6
				Mod.3	3.2	8		82	8.2
				Mod.4	2.4	8		76	7.6
			!Surface	1.4m		[Surface		10m
	f p =20Hz	-	-	i L 7.6m Quarry	Modi 10m		76°	7.6m	10m
								Quarry	Mod.2
			Surface	3 2m		(Surface		2.4m
					10m				10m
				7.6m
							76°7.6m
				Quarry	Mod3
								Quarry	Mod.4

Table 6 :

 6 Different models simulated dynamically.

Table 8 :

 8 Results in dynamic mode with structure.

Table 9 :

 9 Results and observations for the profile of the sinkhole.

		Mode	(m)	Shape of bell h(m) a(°) D(m)
	Without	Static	0.0	6	76	7.6
	structure	Dynamic	2.0	8	78	7.6
	With	Static	1.4	8	76	7.6
	structure	Dynamic	2.5	8	78	7.6
		Observations	2.5	8	80	7.5


Baconnier Simon
email: simon.baconnier-etudiant@ineris.fr

Sidney B Lang
email: lang@bgumail.bgu.ac.il

Seze De

Rene

New Crystal in the Pineal Gland: Characterization and Potential Role in Electromechano-Transduction

The pineal gland is a neuroendocrine transducer secreting melatonin, responsible for the physiological circadian rhythm control. A new form of biomineralization has been studied in the human pineal gland. It consists of small crystals that are less than 20 µm in length. These crystals could be responsible for an electromechanical biological transduction mechanism in the pineal gland due to their structure and piezoelectric properties. Using scanning electron microscopy (SEM) and energy dispersive spectroscopy (EDS), we identified crystals morphology and showed that they only contain calcium, carbon and oxygen elements. Furthermore, the selected-area electron diffraction (SAED) and near-infrared Raman spectroscopy established that the crystals are calcite. We will now focus on the physiological effect of microcrystals in pinealocyte cell culture under Radio-Frequency Electromagnetic-Fields (RF-EMF).

INTRODUCTION

Because of the fast development of mobile telecommunication, the interaction of Electromagnetic Fields (EMF) with biological environment becomes a public health concern. Although the action of non-ionizing radiation on biology is still unclear, several hypotheses of interaction have been suggested: hot spot phenomena, ADN/RF-EMF interaction, EMF effect on cellular development (oncology) [START_REF] Adey | Bioeffects of mobile communications fields; possible mechanisms for cumulative dose[END_REF][START_REF] Maes | 954 MHz microwaves enhance the mutagenic properties of mytomycinC[END_REF][START_REF] Stagg | DNA Synthesis and cell proliferation in C6 glioma and primary glial cells exposed to a 836.55MHz modulated radiofrequency field[END_REF]. But no convincing study brings to the conclusion of an effective risk of RF-EMF for health. The pineal gland converts a neural signal into an endocrine output. The most important hormone it secretes is melatonin the main role of which is to control the physiological circadian rhythm [START_REF] Karasek | Melatonin in humans -where we are 40 years after its discovery[END_REF]. Two biomineralization forms can be observed in the pineal gland. Concretions so called "brain sand", a polycrystalline complex of few millimeters long, and microcrystals the length of which does not exceed 20 micrometers. While concretions have been extensively studied [START_REF] Bocchi | Physical, chemical and mineralogical characterization of carbonate-hydroxyapatite concretions of the human pineal gland[END_REF][START_REF] Galliani | Histochemichal and ultrastructural study of the human pineal gland in the course of aging[END_REF][START_REF] Japha | Calcified inclusions in the superficial pineal gland of the Mongolian gerbil, Meriones unguiculatus[END_REF][START_REF] Kodaka | Scanning electron microscopy and electron probe microanalysis studies of human pineal concretions[END_REF][START_REF] Krstic | Ultrastructural and X-ray microprobe comparison of gerbil and human pineal acervuli[END_REF] no study has been published on the microcrystals. In this article the microcrystals were analyzed with different biophysical techniques. Their physicochemical properties and particularly piezoelectricity would give them an active role in a potential mechanism of electromechanotransduction in the pineal body. We are currently planning a study on the effects of Global System for Mobile (GSM) waves on these microcrystals in cellular culture and their influence on the pineal body physiology.

MATERIALS AND METHODS

The microcrystals were isolated from the pineal bodies using a procedure developed by Weiner and Price [START_REF] Weiner | Disaggregation of bone into crystals[END_REF]. Small pieces of the pineal body (about 10 mg) were placed in a micro-centrifuge tube containing 1.5 ml of 2.5% sodium hypochlorite (diluted commercial bleach) and sonicated for 20 minutes. After allowing the sample to settle for 1 minute, the supernatant liquid was transferred to a second micro-centrifuge tube and centrifuged at approximately 9000 g for 1 minute. The pellet containing the solids was immediately washed twice with 95% ethanol and then resuspended in approximately 50 µl of 100% ethanol. It should be emphasized that, at no point, did any of the samples come into contact with solutions containing calcium ions. SEM samples were collected on transmission electron microscopy grids and analyzed with a JEOL JSM 5600 SEM. Microanalysis studies were performed with a NORAN EDS Analyzing System. Because the microcrystals were initially too thick for High Resolution Transmission Electron Microscopy (HRTEM) observation, they were first crushed between two glass slides. They were then studied with a JEOL-2010 transmission electron microscope equipped with an analytical ISIS system for energy dispersive X-ray spectroscopy (EDS). Near infrared Raman spectra of isolated crystals and pure calcite were obtained with a Bruker IFS 66 FTIR spectrometer equipped with an FRA 106 Raman module and a Ramanscope microscope. Measurements were made with a 40X objective (spot size ~25 µm). The spectral resolution was 2 cm -1 . The samples were excited at 1064 nm using a diode-pumped Nd:YAG laser at about 5 mW power. Second Harmonic Generation (SHG) studies were made with a Nd-YAG laser which produced radiation at 1064 nm and the detected SHG was at 532 nm.

RESULTS AND DISCUSSION

The SEM studies of single microcrystals permitted high-quality morphological analysis. The most common morphology was a very rough cylindrical body with sharp extremities (Fig. 1) that comprised about 95% of the samples observed. The crystal size varied from 1 to about 20 µm. The EDS analyzer coupled to the SEM identified calcium, carbon and oxygen to be the principal elements. Among biominerals containing those atoms, only calcium carbonate and calcium oxalate are potential candidates. The electron diffraction patterns taken from the particles were indexed in terms of a hexagonal unit cell. Near IR Raman spectra were measured on both the microcrystals and on pure calcite powder. The agreement of the peaks was excellent (Fig. 2), confirming the identification of the crystals as calcite (calcium carbonate). We were unable to detect SHG neither in pure hydroxyapatite powder nor in the large pineal concretions. The similarity of the intensity of the SHG in pure calcite to that observed in earlier work on pineal tissue samples [START_REF] Lang | Piezoelectricity in the human pineal gland[END_REF] and the absence of SHG in the large concretions let us think the calcite microcrystals would be the source of the SHG in the previous observation.

The pineal microcrystals appear as a stack of thin rhombohedrons with their flat faces normal to the long axis of the crystal (Fig. 3). These complex structures can be classified using the texture point group nomenclature of Shubnikov et al. [START_REF] Shubnikov | Etudes des textures piezoelectriques[END_REF]. The texture may be noncentrosymmetric because of the structural organization of the sub-unit, even though the single crystals do have a center of symmetry. This symmetry breaking would allow both SHG and piezoelectricity. Calcite in otoconia, microcrystals found in the inner-ear otolith, has been shown to exhibit piezoelectricity [START_REF] Ross | Some properties of otoconia[END_REF][START_REF] Morris | Piezoelectric property of otoliths[END_REF]. These crystals have a structure similar to that of the pineal microcrystals. By that very fact the piezoelectric property of the crystals would allow them to interact with the electrical component of electromagnetic fields. A simplified formula applied to those crystals (f = v/2d) lets us think that these crystals could be sensitive to RF-EMF in the range of 500MHz to 2.5GHz depending on there size. This range contains portable wireless frequencies, GSM (872-960MHz), DCS (1710-1875MHz), UMTS (1900-1920MHz, 2010-2025MHz), or BlueTooth (2400-2483,5MHz). Piezoelectric determination of minute grain requires developing new methods based on either MEMS Precision Instruments microtweezers or direct correlation between electro-optic and piezoelectric properties in crystal with optical microscopy. We introduce a novel approach of the biophysical effects of weak microwave radiation.

CONCLUSION AND PERSPECTIVES

We report here the presence of a new form of mineral deposits in the pineal gland. The calcite microcrystals would have piezoelectric properties with excitability in the frequency range of mobile communications. Their interaction with GSM waves could constitute a new mechanism of electromecano-transduction on the pinealocyte membrane, influencing by the fact the melatonin production. The RF-EMF electrical component interaction with the crystals could induce a morphological modification of the crystals, a vibration depending on the EMF frequency. This morphological change, even tiny, could involve a modification of their cellular environment, by a localized modification of the cellular membrane of related cells. The membrane changes could alter the adrenergic suggested and/or calcium channel function. A similar mechanism of magneto-transduction was revealed by Kirschvink in connection with magnetite crystals of the brain and their interaction with the magnetic component of RF-EMF [START_REF] Kirschvink | Microwave absorption by magnetite: a possible mechanism for coupling nonthermal levels of radiation to biological systems[END_REF]. Pinealocyte can "communicate" through their gap junction [START_REF] Moller | The ultrastructure of the human fetal pineal gland. I. Cell types and blood vessels[END_REF][START_REF] Saez | Pinealocytes in rats: connexin identificatioon and increase in coupling caused by norepinephrine[END_REF]. The deformation caused by the crystal vibrations could thus by simple activation of one or two pinealocytes, activate a whole area of pineal cells and thereby act on the pineal physiology.

The scientific project to be developed is to determine the influence of RF-EMF GSM on pinealocyte and pineal gland physiology through the electromechano-transduction produce by the pineal microcrystals. Using ELISA tests and Laser Scanning Confocal Microscopy we are going to study the evolution in melatonin production and variation in cell calcium flux in primary pineal cell culture.

Fig. 1 :

 1 FIGURES

Fig. 2 :

 2 Fig. 2: Raman spectra of pineal gland microcrystal and pure calcite powder.

Fig. 3 :

 3 Fig. 3: SEM picture of isolated microcrystal on a Formvar-covered TEM grid. Notice the multilayered structure and sharpen extremities. Sketch of two ( :m)T type texture, with or without rotation.


Lynda Priad
email: lynda.driad@mines.u-nancy.fr

Bodgan Piwakowski
email: bogdan.piwakowski@ec-lille.fr.

DETECTION AND CHARACTERISATION OF UNDERGROUND CAVITIES ÜSING HIGH RESOLUTION SEISMIC REFLECTION (HRSR)

INTRODÜCTION

Old Underground mining works are subjected to risks of collapse of which the assessment is a complex process. In the case of inaccessible Underground mines, the difficult task for the risk assessment is to classify the already collapsed zones (no further risk) from still stable structures which may collapse in the future. The knowledge of the characteristics of the Underground cavities and hosting rock mass can provide relevant information. However, accurate Underground cavities detection and characterisation based on geophysical techniques are still a scientific challenge in the subsoil prospecting domain. Among these techniques, the high resolution seismic reflection is the most successfui in the frame of Underground cavities detection. The study presented in this paper, is a part of an extended research program aiming at prediction and controlling ground motions induced by Underground cavities. The N ôbjective is to identify the most \ : appropriate geophysical techniqüe • to L-, : localize and characterize the Underground > cavities at variable depths (several meters to 300 m). Thus, to meet these goals we have investigated the high resolution seismic reflection techniqüe (HRSR). The selected test site is located dose to Greasque municipality in south France (figure 1). In this region, the underground IMI.I K > -r -' "• • ' l -|Xl1 works undertaken during the last two • ) centuries consist of two mines: coal mine \_J and cement stone quarries' (marl limestone). The geological setting is characterized mostly by coal layers and massive limestone formations of the upper cretaceous with altemated coal strata of variable thickness dipping westwards.

V C I Ml.'Ill fIii ii i it

In this paper, shallow mine is assimilated to an Underground quarry

HRSR EXPERIMENT

The HRSR data acquisition 2 was carried out using a Geometrics Strata View Seismograph recording Station, 100 Hz geophones and a silenced "Betsy gun" [l] firing 12-gauge hunting cartridges source. Three profiles A, B and C were carried out at ground level (Figure 2). The line C, purposely positioned mostly on an unmined part, is considered äs a reference and allows the imaging of the natural geological structure. This information is used to interpret the cavities induced anomalies without ambiguities. In addition, this profile crosses two separate drifts in order to check the obtained lateral resolution. The profile B crosses both mined and unmined zones; finally the line A crosses the total area of the quarry in N-S direction, but its CMP (Common Mid Point) resolution is two times smaller which implies worse quality of the seismic image. The data analysis based on the common midpoint method (CMP), provided HRSR 2D sections imaging the crossed Underground structures with a good signal to noise ratio. The geometry acquisition design parameters which have been chosen, regarding specific objectives and expected resolutions are summarized in table l.

HRSR DATA ANALYSIS AND INTERPRETATION

The detection of Underground cavities using seismic sections is based on the analysis of induced signal disturbances. The wave propagation through opened cavities implies time delay and masking effect on the underlying reflectors [START_REF] Piwakowski | Underground tunnels detection and location by high resolution seismic reflection[END_REF], [START_REF] Piwakowski | Imagerie du proche sous sol par sismique haute resolution: Application ä la detection des cavites[END_REF]. In the case where the cavities are collapsed, the seismic response evidences total loss of signal coherence and an important drop in the signal level.

On the three sections of the profiles A, B, C shown on the figure 1, we clearly observe two main features: seismic responses of Underground quarry and unmined structures. The HRSR image of the unmined structures appears äs a succession of cement stone, limestone and coal layers (UM zones). The mined zone is characterized by non observed deeper reflectors (masking effect) and disturbances on the cement stone top layer reflections (MQ zones). In addition, all geological faults met during the cement stone mining are observed. The major faults detected in the northern and southern parts of the quarry (faults Fl and F4) probably implied the limitation of the Underground works in those directions. The profile C evidences unmined zone represented by the geological interfaces of cement stone, coal and limestone layers (CS, C, Lim); in the middle part, a fault (F3) is observed which has been probably the reason of the existing protection pillar (UCP). The seismic proflies A and B offer Images of several cavities effects, in some parts of the sections the almost total loss of the signal suggests that the cavities are collapsed (CC). The roof of the cavities is located at the estimated depth of 15 m. The anomaly CZ is identified on both profiles A and B and is interpreted äs a collapsed zone. This interpretation has been validated and confirmed by the surface observations of the collapse (Rock cracks, ground subsidence). Another interesting seismic feature pointed out on the prbfile A is the relevant consistency with the Underground quarry map, of the detected 15 m width pillar (UCP) structure which appears in the seismic section äs the narrow zone of the continuous and undisturbed reflectors. This observation can be considered äs resolution check oftheHRSR imaging.

CONCLUSION

The high resolution seismic reflection applied to the Underground cavities detection has been tested on the Greasque site (south France), characterised by old coal mines and cement stone quarries. The HRSR enabled imaging the Underground structures at the depth ränge of 7-20 m with quite good resolution. The analysis of the seismic sections evidenced different features and anomalies associated with cavities effect and geological structures. In the quarry zone, the collapse observed at the surface has been clearly identified on the HRSR sections.

In addition, the seismic response of pillars and geological faults correlate quite well with their location on the (rough) quarry Underground map.

Results from this experiment show that the HRSR technique is an interesting and efficient mean capable to evidence the presence of cavities and to recognise their state (already collapsed or not). However, to overcome the experimental limits, further modelling work and borehole check could be undertaken tq refine the characteristics of the cavities and validate the HRSR interpretations.

. ", The experimental feedback of the HRSR investigations in the frame of Underground cavities detection and characterisation confirms its important contribution for subsidence risk assessment. Further investigations are being carried out by INERIS through its research programs aiming at developing an integrated procedure based on geological, mining and geophysical fundamental input in order to develop an elaborated a risk management methodology.

North

Figure l :

 l Figure l: Location map of the test site, Greasque, south France

Table l :

 l Geometry parameters of the seismic profiles carried out in the cement stone quarry

		Length	Geophone	CMP	Exploitable	Lateral	Depth
		(m)	spacing	spacing	frequency	resolution	resolution
			(m)	(m)	band	at 7 m	(V4)
		___			(Hz)	(m)	(m)
	Profil A	335	2	.1	80-180		5.8
	Profil B	105	1	-, 0.5	250-400	3.15	2.3
	ProfilC	60	1 .	0.5	250-400	3.15	2.3

' Conducted by SOVEP, 59338 Tourcoing, France, 52 nie de Roubaix, http://www.sovep.com


Salvi Olivier
email: olivier.salvi@ineris.fr

Gaston Didier

Risk assessment in decision making related to land-use planning (LUP) as required by the Seveso II directive

1 LUP as part of a risk management policy This chapter describe the requirement of the Seveso II directive and an analysis of the use of LUP as part of a risk management policy.

LUP in the Seveso II directive

The article 12 of the Seveso II directive requires that : "Member States shall ensure that their land-use and/or other relevant policies and the procedures for implementing those policies take account of the need, in the long term, to maintain appropriate distances between establishments covered by this Directive and residential areas, areas of public use and areas of particular natural sensitivity or interest, and, Council Directive 96/82/EC of 9 December 1996 on the control of major-accident hazards involving dangerous substances

Introduction

The 1999 annual report from the European Environment Agency indicated that in spite of measures on major industrial accidents in force since 1984, the trend in accidents shows that many of the often seemingly trivial 'lessons learned' from accidents have not yet been sufficiently evaluated and implemented in industry's practices and standards. It was confirmed by the recent accidents (Enschede in 2000, Toulouse in 2001) that have shown that disasters continue to occur throughout the EU despite of the efforts to control major accident hazards with Seveso I and Seveso II directives. Moreover, a recent communication from the Commission (European Commission, 2002) indicates that the frequency of major accidents is estimated at about 3.10~3 per year if we consider the number of accidents reported in the accident database of the European Commission MARS (Major Accident Reporting System) versus the number of hazardous installations throughout the European Union.

So, controlling major accident hazards by reducing the risk on-site is not sufficient to promote a sustainable development for both industry and urban areas in the next decades. It is necessary to organise the settlement of industrial and urban areas with land-use planning (LUP).

This paper presents an approach to implement LUP in the context of the Seveso II directive 1 . In the first part, the requirements of the directive regarding LUP are given, and the use of LUP in a risk management policy is analysed. Then, in the second part, the authors describe the current situation in France, where LUP is applied for more than 10 years with an approach based on the evaluation of the consequences of major accident scenarios. In the third part, the authors share their thoughts to improve the LUP procedure by distinguishing a risk assessment phase and a negotiation phase involving the stakeholders.

in the case of existing establishments, of the need for additional technical measures in accordance with Article 5 so as not to increase the risks to people." Even if the directive is the same within the European Union, the history of the legal systems in the different Member States shows that LUP can be implemented in several ways, using a deterministic approach, a probabilistic approach or tables of appropriate distances [START_REF] Christou | Guidance on land-use planning as required by the council directive 96/82/EC[END_REF].

Risk management policy

To apply the Seveso II directive requirements and to develop a system on the control of major-accident hazards involving dangerous substances, the Member States need to built a risk management policy combining several tools. To describe the risk management policy implemented in France, it is proposed to use the accident model constituted by the system :

-hazardous source : it is constituted by the hazardous installations, products, equipment, processes, -flux of hazard : it is constituted by the dangerous phenomenon like dispersion, fire or explosion, -receptors : they are the elements located in the vicinity of the hazardous installation that could be affected by the flux of hazard. This model is presented in Figure 1.

Hazardous installation

Products Equipment Processes

Hazard flux

Propagation

Receptors

Human Environment Infrastructures

Land-use planning

Hazard potential reduction Risk control

Information of the public

Emergency plans

Figure 1 : Risk management policy in France

The Figure 1 presents also the four components of the risk management policy :

• risk reduction on-site : to reduce the hazard potential (amount of hazardous substances, process pressures and temperature...) and to implement risk control measures during plant operation,

• land-use planning : to maintain appropriate distances between the hazardous installation and the sensitive receptors,

• emergency plans : to prepare the public and the rescue services to behave efficiently in case of a major accident,

• information to the public : to inform the public about the possible consequences of a major accident and the behaviour they need to have in such a situation, but also to explain about the risk reduction on-site and the land-use planning objectives.

Using the model of Figure 1, it is clear that LUP is a complementary tool to the risk reduction on-site and to the emergency plans.

The philosophy behind LUP is to maintain a distance between the hazardous source and the receptors in case of a major accident, in particular if the accident can occur so quickly that the emergencies plans cannot be activated, or so that the number of people to evacuate is too high.

Another important issue is to be aware that LUP do not replace the risk reduction on-site, neither the emergency plans : all these tools are complementary. Moreover, LUP should not be applied before the hazard reduction phase, which is the start of any risk management procedure.

LUP : a decision making process based on risk assessment

As described in the Draft ISO Guide 73 (2001) on risk management, it is fundamental to distinguish the risk assessment process and the decision making process.

The guide explains that risk assessment is a part of the risk management process, ended up with the decision making. Risk assessment is a tool used to measure the risk, characterised by the likelihood and severity of specific events. Risk based decision making process is naturally based on the risk assessment criteria, but must integrate also other criteria that can be cultural, economical, ethical... The decisions related to risk management have to be made for the interest of the Civil Society. The main goal of risk based decision making is to protect people and the environment.

Land-use planning is naturally a decision making process based on risk assessment. Having said that, the recommendations of the ISO Guide 73 should apply to the LUP procedure.

Current situation in France

The risk acceptance of hazardous establishments is a very complex decision process based on several criteria difficult'to evaluate. In France, the only explicit criteria are those related to the consequences of accidents that are used to define the safety distances around hazardous establishments. They are linked to the LUP process. This chapter present some key elements of the history of legal system and the current situation in France.

Obligation to keep the population away from the hazards

In fact, in the Code related to the control of hazardous establishments, the deliverance of the license to operate can be subordinated to a sufficient distance between the establishment and people located around (Code de l'Environnement, Livre V, Art. L512). In other words, the regulatory bodies can not theoretically license new establishments that can potentially harm people in case of a major accident. This requirement is easier to apply for new installations than for existing ones. This idea appeared first in an imperial act in 1810, then in a law dated 1917 on hazardous and unsanitary plants, but it was clearly reinforced with the law n°76-633 dated 19/07/1976 that was integrated in the 'Environment Code' on 18 th September 2000.

Development of the French deterministic approach for risk assessment

To enhance a harmonised implementation of the philosophy of the law presented above, the French Ministry of Environment published in 1990 a guide (Secretary of State to the Prime Minister for the Environment and the Prevention of major technological and natural risks, 1990) that defined reference accident scenarios to be considered for determination of safety distances for land-use planning.

Then, in a circular letter published in June 1992, it was explained that the safety distances should be determined using the safety reports written for the licensing procedure. Consequently, to write the safety reports that aim is characterising the hazard of their plants, the operators have been asked to focus on the examination of the reference scenarios described in the guide. This approach is known as the consequence-based approach or as the French deterministic approach.

The term deterministic means that the assessed scenarios are pre-defined (See Table 1). They are considered independently of their likelihood, which is not assessed. The underlying philosophy is based on the idea that if sufficient measures exist to protect the population from the worst accidents, sufficient protection will also be available for any less serious incident.

Criteria for the evaluation of the consequences of reference scenarios

The criteria used to characterise the hazards are expressed only in terms of effect thresholds. They are used to evaluate the safety distances around the establishments:

-the distance at which the first death occur corresponding to probability of 1 % fatality, -the distance at which irreversible health effects occur.

The effects that are calculated with mathematical models are different according to the physical phenomenon involved in the accident scenario. The effect thresholds correspond to the phenomenon. A list of the thresholds is given in Table 2.

Table 2 : Thresholds for the evaluation of the consequences of the reference accident scenarios

Effects studied

Thermal radiation

Overpressure

Toxic dose

Criteria corresponding to first deaths 5 kW/m 2 (if the exposure is more than 1 min) or heat load of 1000 ( * } (kW/m 2 ) m .s in case of a short exposure duration 140 mbar Based on hC\% and exposure time (passage of the cloud) Criteria corresponding to first irreversible effects 3 kW/m (if the exposure is more than 1 min) or heat load of 600 ( * ) (kW/m 2 ) 4/3 .s in case of a short exposure duration 50 mbar Based on irreversible effects (first injuries) and exposure time (passage of the cloud) (*). : these values are commonly used, but still under discussion It can be noticed that, in case of fragment projection after an explosion, the maximal distance is evaluated but it is often not taken into account for the determination of the safety distance. Harmonised thresholds help to a coherent evaluation of the reference scenario consequences at a national level. Therefore, the French Ministry for Environment has published in 1998 a list of the most frequently used toxic substances (Service de l'Environnement Industriel, 1998). In this list, thresholds corresponding to the first irreversible effects are given. The ongoing work on the determination of the toxic threshold are periodically published on http://www.ineris.fr.

The guidance for the selection of the reference scenarios and the explicit definition of the criteria to evaluate the consequences has contribute to have a consistent procedure, but at the moment we can still notice some discrepancies in the results for two similar installations. Indeed a lot of parameters not defined in the guide have a great influence on the evaluation of the reference scenario consequences, in particular for toxic substances. For example, some parameters should be more precisely defined in the guidance : -The location of the line rupture or leak as a function of the pipe schedule, -The duration of the leak according to safety loop with quick-cut valves, -The duration of the exposure for long leakage, regarding possible response time for people evacuation or confinement, -Meteorological conditions. Moreover, we can still notice some discrepancies related to the models used for the consequence calculation and the assumptions made by the risk analysts.

3 How to improve the situation Aware of the difficulties to implement land-use planning based on the risk assessment carried out for the safety report, the French Competent Authorities in charge of the Seveso II directive have launched a work on this topic in 2000. INERIS has taken part in this work, on one hand, throughout a project entitled "Risk analysis and major accident prevention", and on the other hand, throughout the critical analyses of the safety reports asked by the authorities2 . Thanks to this work, lessons have been learnt about the implementation of LUP and the risk assessment used for the LUP purposes. They are presented hereafter as guiding principles that could be used to improve the LUP procedure in France. However, some of them are so general, that they could be applied in other Member States of the European Union. The views expressed in this paper are entirely the authors. They were partly presented in the report written for the mission of the General Inspectorate of the Environment after the Toulouse accident last year [START_REF] Barthélémy François | Usine de la société Grande Paroisse à Toulouse -Accident du 21 septembre 2001[END_REF], see annex T in tome lib).

Distinguish risk assessment for giving the permit to operate and for LUP

In the safety report required by the Seveso II directive, the risk assessment should enable the operator to conduct a deep analysis on the hazards of a plant and the way to control the associated risks. The analysis should start with a systematic hazard identification and a risk analysis. Then, after the assessment of a set of scenarios, the adequacy of the risk control measures should be discussed. This requires the identification of worst case scenarios and of scenarios that take into account the safety barriers implemented by the operator. The comparison of these two types of scenarios shows the risk reduction and the gain provided by the safety control system. In other words, the demonstration required by article 9 of the Seveso II directive, in particular that "major-accident hazards have been identified and that the necessary measures have been taken to prevent such accidents and to limit their consequences for man and the environment" need a detailed risk assessment. Based on this risk assessment, the authority can verify that the operator has "taken all measures necessary to prevent major accidents and to limit their consequences for man and the environment", as required by the Article 5 of the Seveso II directive.

The purpose of LUP is, as explain in Article 12 of the directive to "maintain appropriate distances between establishments [...] and residential areas, areas of public use and areas of particular natural sensitivity or interest". This objective can be interpreted as the need to estimate the consequences of major accidents in order to determine the zones in which the number of people should be minimized. Having in mind this objective, it is proposed to estimate the consequences of accidents to define the zoning around hazardous establishments. This approach allows to reason on the hazard potential and not on the risk (combination of the severity and probability of the accidents). It is important for 2 reasons :

-risk control measures can fail, then it is useful to know to what extend damages can occurred in case of an accident, there is less uncertainty in consequence assessment than in risk assessment. To assess the consequences of accidents, scenarios must be defined. Since an accident scenario is the combination and the conjunction of several events, the assumptions related to the conditions of occurrence of the events including the functioning or not of some safety devices have a great importance on the final results.

Scenarios taking into account safety barriers

The scenarios that can lead to the most serious consequences are usually used to determine the emergency plans, The scenarios considering that the safety barriers in place can prevent the occurrence of some events or reduce their consequences can also be identified. They are usually used to judge the adequacy of the safety control system which includes naturally the safety equipment and operations, and the organisational measures, both defined and animated by the safety management system. For land-use planning purpose, intermediate scenarios are needed. Such scenarios should take into account some safety barriers, for example, those prescribed by standards and regulations. In particular, the safety barriers that are the most reliable should be considered. These safety barriers should be identified, reviewed and consensually accepted by the authorities and the other stakeholders of the LUP process.

Defining certain assumptions related to the consequence assessment of the scenarios

In France, considering several risk assessments of similar hazardous installations, strong discrepancies appear on the estimation of the distances related to the apparition of the lethal effects and first injuries. This statement was the same in the ASSURANCE project [START_REF] Lauridsen | Uncertainties in risk analysis of chemical establishments -the ASSURANCE project -1998-2001[END_REF]Makert et al., 2001). This benchmark exercise allows to measure the size of the uncertainties in risk assessment and to propose some explanations. The results indicate discrepancies of about 4 orders of magnitude (variation of 1.5 10" 7 to 2.1 10 3 ) in the probability estimation of some events. There is a factor 6 for the estimation of the consequences of a liquefied ammonia release in unfavourable atmospheric conditions (variation from 1510 m to 9700 m for first injuries, and from 570 m to 3800 m for the first deaths). These results are due to the great variety of possible assumptions related to the parameters describing the conditions of occurrence of the scenarios and the hazardous phenomena. The choice of these assumptions generates uncertainties on the estimated safety distances.

However, the decision making process constituted by the LUP should be homogeneous and consistent at the national level, and further at the European level. Therefore, defining certain assumptions could reduce the discrepancy in the decision process, and should make it more transparent. It is no more tolerable that the results depends on the experience, the feeling of a risk evaluator, neither the pressures related to the implication of the results of its estimations. It is essential that the rules are the same for all, should they be discussed and determined after a debate involving all stakeholders : industry, authorities, mayors, communities and risk appraisers.

In order to make it explicit that the scenarios used to estimate the danger zones in LUP are based on predefined assumptions, it is proposed to call them conventional or consensual scenarios.

LUP, local concerns

LUP is implemented at the local level and should take account of local specificity. The LUP ACS project [START_REF] Duijm | Land Use Planning And Chemical Sites -the LUPACS project -1996-2000[END_REF] proposed a methodology for land-use planning involving chemical sites for making decision in local and regional administration. This European project expressed that LUP is a complex process involving actors at different decision making levels with different interests. In that context, the estimation of the consequences of the proposed consensual scenarios should be one of the criteria used to define the danger zones for the use of land. But the problem has other dimensions than the risk level, and decision makers should take into account the social, economical, cultural dimensions of the problem. The danger zones based on the consensual scenarios enable to identify the residential areas and other sensitive areas or infrastructures capable to be affected by the accident. They identify also the stakeholders who should be involved in the discussions and negotitation that will lead to find a consensus on the most adequate alternatives between :

• keep large zones around the plant and disable the development of the city;

• reduce the zones by implementing technical measures to enable the development of the city.

Both for new installations or existing ones, the danger zones should not be directly transposed. The implication of land-use by taking into account the hazard potential of a plant should be discussed at the local level, and should be analysed regarding all the aspects of the decision.

Suggestion of a 2 steps approach

Using the lessons learnt for more than 10 years, INERIS promotes a two-step approach that allows :

• To assess the hazard and estimate the danger zones in a homogeneous manner, following a commun procedure ;

• To develop and assess alternatives that can be negotiated at the local level, according to criteria that should be made explicit.

This approach is presented in Figure 2. The first part of the approach is the risk assessment that leads to the estimation of the danger zones and their communication to the public, in order the public and the stakeholders know about the extension of the damages in case of a severe accident. This step is based on the examination of a set of scenarios in the safety report, covering the different hazardous phenomenon associated to accidents involving toxic, explosive and flammable substances. The consensus scenarios used to determine the danger zones should be defined in a guide, which is currently under development.

The second part is the treatment of the risk and the negotiation of the alternatives for the landuse planning, according to the local constraints like the composition of the environment in the vicinity of the hazardous installations, and socio-economical context. At this second stage, the discussions and the negotiation should associate the various stakeholders identified within the danger zones determined in the first stage. This phase of the procedure can be iterative if the local constraints need the implementation of several risk reducing measures to reduce the zones. However, attention should be paid in order to avoid that the final danger zones become too small, because then the principles of LUP (maintain a appropriate distance between the hazardous installations and the receptors in case of a major accident) would not be respected. This two-step approach bring consistency in the hazard assessment, and transparency in the negotiation so as flexibility to take into account the local context.

Conclusion

As written in 1999, in the European Environment Agency report, major accidents continue to occur despite of the efforts of the industry and of the authorities. "Zero risk doesn't exist" is a reality. This motto should reinforce the conviction that all stakeholders, in particular, the industry, the authorities and also the civil society should optimise the control of the risk in a best way, as a function of the resources of each, and of the benefits associated to the industrial activity.

It is also important to underline that the choice of accident scenarios for LUP purposes doesn't mean that a more catastrophic accident will not happen. The choice of scenarios and the estimation of their consequences is a tool to make a political decision.

The two-step approach described in this paper results from the experience of more than 10 years of application of LUP in France. It is based on an assessment of the consequences of accident scenarios and on transparent negotiation of the risk reducing measures to determine the final danger zones that are inscribed on the urbanism documents. This suggested approach brings some elements of response to the European Council that expressed, at its Barcelona meeting, the invitation " to develop a strategic approach on the management of technological risks, while considering the social, economic and environmental challenges in relation with the issue of sustainable urban development."

  Figure 2 : Two-step approach for land-use planning

Table 1 :

 1 Reference accident scenarios to be considered for land-use planning purposes

	Type of risks and facilities	Type of accident scenario
	Risks linked to liquefied combustible gas	Scenario A :
	facilities (fixed, semi-mobile or mobile)	BLEVE (Boiling Liquid Expanding
		Vapor Explosion)
		Scenario B :
		VCE (Vapour Cloud Explosion)
	Risks linked to vessels containing liquified or	Scenario C :
	non-liquified toxic gases where the containment is not designed to resist external	Total instantaneous loss of containment
	damage or internal reactions of products
	Risks linked to vessels containing toxic gases	Scenario D :
	where the containment is designed to resist external damage or internal reactions of products	Instantaneous rupture of the largest pipeline leading to the highest mass flow

The French legal system allows the Authorities to ask a third party expert to make a critical analysis of a safety report provided by an operator of a Seveso plant. The critical analysis is paid by the operator and is delivered to both him and the Authorities.

ACKNOWLEDGEMENT

This paper was partly supported by the project mentioned in the text that is funded by the French Environment Ministry.