INFLUENCE OF MINING EXCAVATION ON ENERGY REDISTRIBUTION AND ROCKBURST POTENTIAL ABSTRACT

Keywords: Deep mine, Numerical Modeling, Rockburst, Energy balance

With the increase of mining depth, rockburst is one of the most serious disasters which threaten the safety of mine operators and the surface stability. Several approaches have been developed since the sixties to assess rockburst potential in underground hardrock mines. Some of the approaches are based on energy balance around mining excavations such as the Energy Release Rate (ERR) that was developed in South Africa and more recently the strain Energy Storage Rate (ESR). This paper presents the results of a detailed numerical modeling case study for the assessment of rockburst potential in a room-and-pillar copper mine in Poland. The method of numerical simulation is a quasi-static finite difference method (FDM -FLAC3D). The Energy Storage Rate (ESR) is numerically calculated to predict burst occurrence.

INTRODUCTION

The assessment of rockburst occurrence is essential for the design and the safety of a deep mining exploitation. This type of dynamic event corresponds to a sudden release of energy from a highly stressed rock. Most frequently, rockburst causes rocks to crush and collapse together. Several factors are known to have an impact on the rockburst hazard [START_REF] Butra | Rockburst hazard evaluation and prevention methods in Polish copper mines[END_REF]: depth of mining, lithology of rockmass, thickness of deposit, presence of discontinuity planes, geo-mechanical characteristics of rock-mass. Some of the mining factors also have an influence on the occurrence of rockburst: mining method, roof control method, pattern of deposit cut, concentration of mining operations, spatial limits of mining operations. For mining areas, numerical modeling is often used to identify zones of high stress concentration. Areas where calculated stresses are predicted to exceed the strength criterion of the rock can then be monitored to prevent instabilities. Another approach is the consideration of energy instead of stresses. Various methods involving energy balance have been developed; their main advantage is their consideration of rock's stress and strain state in a unique value of energy. These methods rely on the definition of a critical energy threshold corresponding to the amount of energy that can be stored in the rock before being brutally released as energy waves. For example, the evaluation of the Energy Storage Rate (ESR) provides an indication on the state of the rock surrounding an excavation.

After a review of the methods of detection of mining induced seismicity, this paper presents a numerical modeling study using FLAC-3D software, which is used to calculate the ESR. The efficiency of the modeling technique is demonstrated for a cylindrical excavation in an infinite elastic rockmass under hydrostatic stresses. The calculation of energy is then applied to a case study of a Polish underground mine, showing the need to define accurate levels of critical energy. Finally, a method for the calculation of critical amount of energy stored in the particular case of a triaxial test is presented.

ENERGY BALANCE IN MINING

Review of rockburst detection/prevention methods

The prediction of the occurrence of rockburst is based on detection and prevention methods used by mining companies. Primarily, micro-seismic monitoring systems (consisting of uniaxial or triaxial accelerometers) have become an integral part of most hard rock deep mines in an effort to characterize mining induced seismicity [START_REF] Trifu | Use of microseismic monitoring for rockburst management at Vale Inco mines[END_REF][START_REF] Archibald | Development of in-situ rockburst precursor warning systems[END_REF]. The location, timing and magnitude of the induced seismic events can be determined, but the cause of these events is not directly identifiable. In fact, this seismicity can be created by high stresses, fault slip, fracture propagation, or strain burst. The main objective of monitoring is to adapt ground control support or to calibrate numerical modeling.

Optimization of the mining sequence in order to store high stresses into the free faces and good ground support practices have also been applied to reduce the proneness to rockburst.

The occurrence of strain bursts, a type of rockburst caused by the burst of a pillar, a face or a floor in a mine, has also been investigated based on analytical and numerical modeling. Various methods have been developed to locate the areas that are prone to rockburst occurrence and rockburst indicators were introduced depending on stress concentration (Tadjus, 1997), energy concentration [START_REF] Itasca Wang | Comprehensive prediction of rockburst based on analysis of strain energy in rocks[END_REF], ratios of stress/strength [START_REF] Mitri | Numerical studies of destressing mine pillars in highly stressed rock[END_REF].

Among them are the Energy Storage Rate (ESR) [START_REF] Mitri | FE modelling of mining-induced energy release and storage rates[END_REF] and Energy Release Rate (ERR) [START_REF] Cook | Rock mechanics applied to the study of rockbursts[END_REF], which were introduced to relate rockburst occurrence and the amount of energy stored in rockmass. The use of ERR as a measure of underground conditions rather than an indicator of seismicity was discussed. The present study uses ESR to characterise the amount of energy that can be stored in a rockmass.

Description of the different energy components

As first described by Cook (1967), reviewed by [START_REF] Salamon | Energy considerations in rock mechanics: fundamental results[END_REF], and then refined by Brady and Brown (1985) and [START_REF] Hedley | Rockburst Handbook for Ontario Hardrock Mines[END_REF], every transition from one state of equilibrium to another during a mining operation is associated with energy transfers. The energy components associated with a transition between these arbitrary states of mining are described as follows:

is the work done by external and body forces, is the elastic strain energy initially stored in rock mined during the transition between the two states, is the increase in stored energy, is the released energy, is the work applied on the support of the excavation, and is the related to plastic deformation. The energy balance is therefore defined by:

( + ) - + + = > 0 (1)
An incremental approach is used to follow the changes due to mining (see figure 1). The mining sequence is decomposed into different stages of excavation and the value of ESR is calculated at each stage. Hereafter, the subscript k represents the number of the current stage of excavation, k-1, the number of the previous one, and .., .., ! , " .., " .., ! , the corresponding components of the stress and strain tensors in a Cartesian coordinate system:

### = $ , %, &, %, %%, %&, &, %&, &&, ' " ) = $ " , " %, " &, " %, " %%, " %&, " &, " %&, " &&, '
As noticed by [START_REF] Mitri | FE modelling of mining-induced energy release and storage rates[END_REF], an increase in stored energy in a volume V, can be separated into ! , which is the energy due to induced stress and * , which is the energy due to in-situ stress (initial stress before excavation).

, = !, + *, (2)
Where:

!, = ! * * , ∆" ##### : /01, ####### * 23 , *, = , ∆" ##### : ! ###### * 23 (3)
∆" ##### is the transposed incremental strain tensor associated with the stage k

/01, ####### is the induced stress tensor /01, ####### = ### - ! ######
A FISH command was created in FLAC-3D for the purpose of this study, using a formula adapted to the Cartesian coordinates, applied in each zone of the model (4_678 is the volume of the zone considered): The Energy Storage Rate (ESR) is then given by @AB = 2 23 !, + *, and @AB = F @AB G (4)

Figure 1: Energy components at mining step k (adapted from [START_REF] Mitri | FE modelling of mining-induced energy release and storage rates[END_REF])

Cylindrical excavation in elastic rockmass

We first compare the results of a numerical modeling using the commercial codes UDEC and FLAC-3D with the analytical solution calculated by [START_REF] Salamon | Energy considerations in rock mechanics: fundamental results[END_REF] for a cylindrical opening created instantaneously in a linearly elastic material. The results of this modeling ensure that the value given by the FISH-command (developed in FLAC-3D) is accurate.

Stage 1 corresponds to the excavation of a 1m-radius cylinder and stage 2 corresponds to a 2m-radius excavation (figure 2). The analytical solution can be easily found from this geometry. Work done at the outer boundary is calculated at a distance of 10 m from the center.

Bulk modulus (MPa)

38.9*10 I Shear modulus (MPa) 29.7*10 I Initial stress (MPa) 100 The absence of values for ! and * in the column corresponding to UDEC is explained by the fact that this software does not distinguish between stored energy resulting from induced stress versus in-situ stress. The results are very close to those given by the theory. The predicted relations for stage 1 should be:

J = 0.5 * J = 1 + 2 (5)
These relations are observed with an error of less than 4% and imply that half of the boundary work is stored into the rock mass. As described by [START_REF] Hedley | Rockburst Handbook for Ontario Hardrock Mines[END_REF], this proportion is not the same for the enlargement of the cavity. The other part of energy should be dissipated by damping of the spherical waves emitted during the instantaneous excavation.

The same type of calculation was also performed for a spherical excavation in an infinite, elastic and homogenous material, and was compared with the results of the analytical solution and the boundary element program Examine-3D (Rocscience). Again, the results fit the theoretical solution, which allows us to affirm that the FISH command can efficiently predict the amount of energy stored following an excavation sequence while modeling an elastic rockmass in three dimensions.

CASE STUDY

Presentation of the Rudna Copper Mine

The Rudna mine is one of three mines in the Polish Legniga Glogowski Copper Belt, along with the Polkowice mine and the Lubin mine. The mined ore corresponds to a thin and flat rock layer, at an actual average depth of 1250 m below ground surface, exploited since the 1970's with different variant of room-and-pillar method. The production of copper ore in this mine amounts to about 12.5 Mt per year, with a grade below 2%.

Because of the mine depth, the geology and the presence of fault, this mine is subject to many rockbursts. As described by Butra, (1998), even if the geological conditions above this mine are fair (thick and strong dolomite layer overlaying the ore), the size of mined out area (about 75 LM * ), the number of dynamic events and the high stress concentration have made the mining environment increasingly complex. Every year, about 3000 induced tremors are registered in the LGOM district, with the energy of tremors reaching up to 10 !N Joules [START_REF] Orzepowski | Evaluation of rock-mass state in Polish copper mines through monitoring the borehole deformations[END_REF].

Rockbolting and hydraulic backfilling is commonly used to support the roof and to improve the stability during the first stage of excavation. Moreover the method of the yielding remnant pillars is used. This approach provides immediate roof strata protection in the working zone, as well as advantageous stress distribution due to instant support of main roof strata on yield pillars.

Control of this method constitutes a real achievement for the Polish copper mining industry. The main idea behind room-and-pillar mining with a roof-deflected system consists in the creation of pillars whose size permits yielding simultaneously with lowering and deflection on the roof strata. Crushed yield pillars of wedge-column shape behave similarly to a deformable artificial support which does not accumulate elastic strain energy. For the modeling, the pillars were designed to be rectangular instead of wedge.

Geometry, properties and assumptions of the model

According to [START_REF] Katulski | Measurement of pre-mining stressfield at Rudna Mine[END_REF], horizontal principal stresses are considered to be 40% larger than the vertical stress, which is estimated to be about 30 MPa. Model geometry (figure 3 and4) is taken from the map of the Rudna mine provided by KGHM (the company that operates the mine). Attention was focused on a block currently under excavation. The mining method is a one-phase room-and-pillar of type J-3S, the remnant pillars have a width of 8 m.

Geologic strata

Thickness (m)

Description of the mining sequence

The mining sequence followed in the mine is decomposed in 3 stages represented in figure 5. The dimensions of the remnants pillar and their shape shown in the figure 5 Steps A correspond to the development work, required to prepare the flat bedded deposit for mining and also to obtain a precise geologic description of the orebody. These stopes are roadways for ore transport and communication are established inside production stopes. Steps B are the ore recovery steps, corresponding to the creation of the remnant pillars. The backfilling stage by roof caving method is not in the scope of this study.

Results of the modeling

Comparison of the vertical stress distribution with solution given by SUIT 3D

A modeling of the mine was conducted with a boundary element software, SUIT-3D. Similarly to the FLAC-3D result, the distribution of vertical stress shows a concentration of high stresses (> 90 MPa) in the center of mining area, where remnant pillars are created (figure 6 and7). The construction of the geometry with SUIT-3D is easier and the time of calculation is shorter. However, the abilities of FLAC-3D to change the behavior of the material into elastic-plastic with strain hardening/ softening, or to insert faults, dykes and other geological parameters, make FLAC-3D the preferred code for further analysis. Moreover, the calculation of energy components is harder to implement with SUIT-3D. During the first mining steps, the energy stored in the pillars around the excavation increases, especially in the pillars at the border of the initial block. These pillars have already been disturbed during their creation and therefore this perturbation may increase the probability of strainburst occurrence. The creation of the remnant pillars leads to an increase in stored energy. The evolution of the mining sequence shows a "flow" of stored energy following the mining front.

Between 1990 and 2010, the high energy tremors caused 323 events with consequences in the mines of the Legnica-Glogow Mining Basin (Butra, 2011). The number of events could be partially explained by a quantity of energy stored in the pillar close to a critical energy. We describe a way of calculating a critical value of the ESR to be applied in further work that would assess the probability of failure of the pillars. To continue this work, a model with accurate data and a more realistic mining sequence could also be realized.

The analytical solution for the total amount of energy stored in the elastic state at a timestep " < ! " is:

, = mB N * N ( * + (24 -1) =0h * -44 * =0h )
The solution for the total amount of energy dissipated by plastic work in the plastic state at a timestep > ! is: = mB N * 3 ( o -p =0h -2 * J * qo )

Where: o = !rstu (v) ! >0(v)
and p = !rstu (w) ! >0(w)

.

The solution given by the FISH command in FLAC-3D can be verified by studying the theoretical solution in the case of a triaxial test (figure 10). Subsequently, the critical energy density value in the sample can be evaluated for different confining stress (figure 11). The critical energy density is defined by: = , 37 xM (7)

When a zone has its stored energy that reaches this critical energy, it is possible to consider that the pillar is close to failure.

CONCLUSION

This paper focuses on the prediction of mining induced strain bursts using an energy approach. The distinction between the different energy components involved in a mining excavation is discussed. The calculation of these components is implemented into the commercial numerical modeling code FLAC-3D and tested through a comparison with the analytical solution for a cylindrical excavation in an infinite, homogenous and elastic rockmass. It was then applied to an underground copper mine in Poland.

The evolution of the stored elastic energy was considered at each stage of mining. A triaxial test on an elastic-plastic sample was modeled in order to calculate the value of a critical energy density including the influence of the confining stress. The next step of study will be to compare the critical energy to the energy stored in the pillars of the Rudna mine in order to determine which pillars are the most likely to be unstable and to release a high quantity of stored energy as rockbursts.

  * ( >01, , * ∆" , + >01, ,%% * ∆" %%, + >01, ,&& * ∆" &&, + 2 * >01, , % * ∆" %, + >01, , & * ∆" &, + >01, ,%& * ∆" %&, ) *, ,9:;< = 4 &=0 * ( , ! * ∆" , + %%, ! * ∆" %%, + &&, ! * ∆" &&, + 2 * %, ! * ∆" %, + &, ! * ∆" &, + %&, ! * ∆" %&, )

  corresponds to another part of the mine which is out of the scope of the study. The cut A-A shows the first stage of development, the cut B-B the creation of the remnant pillars and the cut C-C a stage of caving of the blocks.

Figure 3 :

 3 Figure 3: Horizontal view and dimensions of the model

Figure 5 :

 5 Figure 5: Examples of stages in the mining sequence

Figure 7 :

 7 Figure 7: Vertical stress at the end of the mining sequence with SUIT-3D (the scale is not the same)

Figure 8 :

 8 Figure 8: Energy storage rate: ESR (PQ/S T ) during a mining sequence in the Rudna Mine

Figure 10 :

 10 Figure 10: Analytical and calculated evolution of stored elastic energy and Plastic work during triaxial test

Table 1 :

 1

			STAGE I		STAGE II
		Analytical	UDEC	FLAC-3D	Analytical	UDEC	FLAC-3D
		value			value
	! (MJ) * (MJ) (MJ)	0.0266 0 0.0266	No value 0.0264 No value 0.002 0.0266 0.0274	0.0581 0.0386 0.0969	No value No value 0.0934	0.0591 0.0385 0.0976
	(MJ) 0.0535	0.0527	0.055	0.157	0.152	0.163

Properties of the material Figure 2: Dimensions of cylindrical excavation model

Table 2 :

 2 Energy balance for a cylindrical excavation in an elastic material

Table 3 :

 3 Properties of the material (adapted from[START_REF] Pytel | Rock mass-mine workings interaction model for Polish copper mine conditions[END_REF]

	Dolomites Orebody	30 3	γ (MN/m3) 2.7e-2 2.7e-2	0.24 4 0.18	E (MPa) 8550 2039	Rc (MPa) 136 67.59	Rt (MPa) 8.8 3
	Sandstones	30	2.7e-2	0.13	15000	18.4	0.8

ACKNOLEDGEMENTS

This study was sponsored by the I 2 Mine European R&D Program and was conducted at INERIS in collaboration with Ecole des Mines de Nancy and McGill Mine Design Laboratory. The authors are grateful for the support provided by the project "Ressources 21" (Resources of the 21 st century) held by the Université de Lorraine. Thanks are also given to KGHM CUPRUM for the data provided.

CALCULATION OF A CRITICAL ENERGY

As reported by [START_REF] Mitri | FE modelling of mining-induced energy release and storage rates[END_REF], [START_REF] Mitri | Assessment of horizontal pillar burst in deep hard rock mines[END_REF], and by [START_REF] Xie | Energy analysis and criteria for structural failure of rocks[END_REF], a value of critical energy has to be determined, representing the maximum amount of elastic energy that can be stored in a rock before it fails. Previous works chose this value based on the uniaxial compressive strength, which is a conservative assumption that does not take into account the effect of a confining pressure. The state of the stresses in the pillars was considered to be uniaxial because the failure is taking place at the pillar skin, where the normal stress must be zero. If we consider the yielding of the center of a pillar, this normal stress is not zero and acts as a confining stress. During a triaxial test, the value of critical energy highly relies on the value of the confining stress. The critical value of energy stored can be obtained by considering that the plastic yielding will lead to an unstable state close to the failure of the pillar.

A triaxial test was modeled with FLAC-3D to determine the critical amount of energy that can be stored depending on the confining pressure. The energy balance during a triaxial test is conducted at a constant axial strain rate (constant downward speed of the upper plate V represented in figure 9). A yield criterion of Mohr-Coulomb type is considered. The material is characterized by its cohesion c, its friction angle U and its dilation angle ψ. The dimensions of the sample are given in figure 9. Let Ω be the volume of the sample and dΩ its surface.When the sample is compressed, a part of the work done at the external boundary is stored while the other part is dissipated by plasticity. The energy balance is given by:

is the boundary loading work = ∭ :

) " ### 23 is the energy dissipated by plastic work


S Coccia

A Lizeur

P Bigarré

E. Klein INERIS I Contrucci

ACCURATE 3D LOCATION OF MINE INDUCED SEISMICITY IN COMPLEX NEAR-FIELD UNDERGROUND CONDITIONS

Installation of local high resolution microseismic networks in the vicinity of hazardous mining fronts, and completing a preexisting mine scale seismic network, is supposed to help detecting much smaller events and processing seismic data with much better accuracy, then gaining insight in the quantification of the rock mass response versus time, space and mine production, all critical parameters for rating the seismic hazard on a continuous basis.

However, the detection and processing of small magnitude events require some precautions due to near-field conditions adverse to the detection of significant changes in the seismic pattern. High resolution monitoring raises new issues related to the magnified complexity of the rock mass when dealing with higher frequency seismic waves travelling through the host rock, intersected by faults and geological disturbances, adjacent backfilled works intertwined with multilevel fast advancing mining works. This magnified complexity may introduce artefacts in the accurate location of the small magnitude sources, blurring the interpretation of the seismicity.

To face such situations of complex mining underground conditions, the authors intend to develop a new approach based on the implementation of 3D dynamic velocity model enabling to take into account not only the geological features surrounding a mining block, but the dynamic mining process itself, i.e. the creation of mining voids and surrounding disturbed zones.

Besides synthetic numerical tests that have been run to assess the relevance of this issue, the research work is tested through the back analysis of an intense microseismic swarm recorded during the brutal caving process of a solution mined cavern. This paper presents this numerical procedure currently being developed for operational implementation in the near future in deep mining works within the strategy of deployments of mobile local acoustic/micro seismic arrays.

INTRODUCTION

Nowadays, the increasing need of ore extraction leads to the use of massive underground mining systems, such as cut-and-fill, sublevel and longhole stoping (artificially supported) or block and panel caving (naturally supported). They are highly productive mining methods, where ground stability issues can be severe and difficult to control for mining operators. For instance, block and panel caving miners have virtually no visual access to the cave and limited opportunities to install geotechnical instrumentation in the working area. In terms of risk mitigation this situation has been improved with the introduction of high resolution seismic system surrounding closely an area of interest (Handley 2006 ).

Since its inception in the 1970s and its fast development from 1990's, seismic monitoring has proved to be invaluable tool for identifying underground areas at risk. However, the traditional mine-scale monitoring with its wide sensor spacing allows neither the recording of small events, nor the precise location of larger ones, limiting strongly the insight of any local seismogenic zone, where the stress gradient ahead of the mining front is breaking the rock (Potvin and Hudyma ).

In addition, in the last decade, the mining production in underground hardrock mines has been also pushed to deeper levels, where the in situ and mine induced stresses are getting higher. This gives emphasis to the increasing necessity to ensure the safety of mine workers and minimize production losses by reinforcing the identification of rockburst potential areas.

To manage this seismic risk different engineering approaches have been proposed (Hudyma and Potvin 2010 ), all based on analysis techniques of high resolution microseismic data.

It is well known that high resolution seismic monitoring can provide great benefits, but improving microseismic monitoring efficiency and accuracy at a mine requires a particular attention and sophisticated studies. With the installation of the microseismic monitoring system, a large amount of data was collected and it provides a unique opportunity for real time characterizations of fracture/failure processes in terms of event location, magnitude and source mechanisms (Iannacchione et al. 2003 ).

Advances in seismic instrumentations have to correspond in advances in techniques of mining-induced seismicity. Without forget that every mine has its own specificity, generating particular type of ground motion, depending on the type of ore, technology and local geology (Zembaty 2004 ).

A good understanding of mine induced seismicity relies strongly on the detection and correct location of all events in order to quantify the rock mass response to both the mining process and the geological mapped structures in the vicinity. Considering the complexity of the rock mass response, seismic hazard may be rated for a mining block as an averaged index over a given time window but still considered to vary rapidly if not brutally, due to heterogeneity of the host rock, presence of highly stressed zones, goaf generation (that part of a mine from which the mineral has been partially or wholly removed; the waste left in old workings; -called also gob), etc., requiring to monitor the most complete induced seismicity catalogue. However, most mine-scale seismic networks are deployed at a much larger scale compared to a single mining work and have then the severe drawback of not detecting the numerous small magnitude microseismic events. This drawback can be managed in some specific situations by deploying locally a high resolution array, completing the coverage and resolution needed for any specific local mining work rated at risk. This may be seen as a temporary array set up and moving along with the mining.

The use of such mobile monitoring systems, aims first to complete largely the seismic catalogue for the area of interest and second to increase drastically the correctness of the calculated position of the high magnitude events, enhancing then significantly the ability to the mining staff to detect any significant change in the seismic activity versus the local geology encountered and the production rate. Concerning the latter, the correct location of any seismic event is critically important because not only a good location provides a guide for a rescue effort and other emergency plan, but all subsequent seismological processing relies on the event position. Eventually, any rockburst mitigation strategy originates from the correct 3D mapping of the seismic sources relatively to all the old or active mining works present in the area at risk. Source location uncertainties too large compared to the different emissive components of the monitored area affect obviously the capabilities to detect those changes ground engineers or mine seismologist are looking for.

It is well known that source location accuracy is affected by many factors. The most important ones are: sensor spatial distribution, accuracy in picking arrival times of seismic phases, knowledge of the velocity model, the computing of ray paths along with the location algorithm. For the mine environment, there are also other practical factors which play a fundamental role, such as mine layout, geology and various ground control issues.

The focus of this paper is to describe a new accurate microseismic event location approach to be used with high resolution monitoring systems in situations in which an evolving 3D velocity model is supposed of importance. This works is being carried out thanks to the I2Mine (Innovative Technologies and Concepts for the Intelligent Deep Mine of the Future) FP7 European project, which focuses on the development of innovative technologies suitable for deep mining activities. It targets more especially health and safety issues, questions related to the working environment and the environmental implications and impacts.

LOCALISATION PARAMETERS

Data used to constrain event locations are usually derived from seismograms recorded at seismic stations, distributed around the monitored area and include arrival times and polarization angles.

Once the data quality is checked and phases are correctly picked the largest influence on the location of seismic events is related to the use of an advanced source location algorithm and the accuracy of the velocity model used for the rock mass. Thus, we have decided to deal with these aspects and their characteristics are discussed in the next sections.

LOCALISATION ALGORITHM

There are two broad groups of event location methods, namely "absolute location" and "relative location" methods. Absolute location is determined or specified within a fixed, geographic coordinate system and a fixed time base (e.g., Coordinated Universal Time, UTC). A relative location uses some or all of the seismograms from two or more neighbouring seismic events to locate them relatively 1980 1992 to each other (Spence and Pavlis ). This location is old and highly popular method, above all because the travel time anomalies resulting from velocity model uncertainties are removed implicitly. However, improvements in relative location accuracy obtained using the double-difference algorithm (Waldhauser and Ellsworth 2000 ' 2002 ) and source specific station terms (Richards-Dinger and Shearer 2000 ) often produce a dramatic sharpening of seismicity patterns, particularly when more accurate timing is obtained using waveform cross correlation (Shearer ). Evaluating the performance of these methods is complicated, by the fact that the true earthquake locations are unknown.

Usually, an earthquake location is determined by an inverse problem: the match or misfit between observed arrival times of seismic wave-energy at seismic stations, and predictions of these arrival times for different source locations using a given elastic-wave speed model (Tarantola and Valette 1982 ).

The general probabilistic solution to the inverse problem of event location from the available data is given by:

(1) ju (d,m) where m denotes the vector of source location parameters (spatial coordinates and origin time) which takes values from parameter space M; p(m)

is the prior pdf (probability density function) representing all information available about the location before (prior to) using the data dobs (e.g., arrival times); /?(d) is the pdf over a data space D describing the data uncertainty in d obs due to measurement and processing uncertainties and ,F(d,m) is the ability of the forward problem (e.g., travel time calculation) to predict the observed data (e.g., arrival times). The constant k normalizes Q to unit integral over DxM and / u(d,m)is the homogeneous distribution over data d and parameters m .

It is usual to call the integral in (1) the likelihood function L(m), which gives a (non-normalized) measure of how good any model m is in explaining the observed data With some simplifications (Lomax et al. 2009 ) a maximum likelihood origin time, t 0 , can be determined analytically from weighted means of the observed arrival times and the predicted travel times (e.g., Tarantola and Valette 1982 ), and if the observed and predicted times are uncorrelated the likelihood function is:

,[yr-?T(x)] 2 (2)
where x is the spatial part of m, T° are observed travel times, Tf are the calculated travel times for observation i (i.e., Tf represents the travel time, rather than arrival time, part of /(m)), and u i summarizes the associated standard deviation of uncertainty in T° and Tf .

Though not normalized, L(x) is sufficient to provide the relative probability of any location m being the best estimate of the event location given the available data measurements. Since in practice integrating over all of DxM to find normalizing constant k in Equation 1 is often computationally intractable, the product of the prior, spatial location information p(x) (i.e., the spatial part of p(m) ) and the non-normalized likelihood L(x) is usually taken as the objective function for inversion and searching in direct search location algorithms. If L(x) is determined throughout the prior pdf p(x) through a global-search, then Equation 1 can be normalized approximately after location. Therefore, according to Lomax et a/. 2005 ' 2009 we refer to such an approximately normalized function, p(x) L(x), as a location pdf.

Moreover, to solve the inverse problem for seismic source location we used the same procedure described by Lomax et al. 2009 and we focused on the Non-Linearized or Direct-Search methods.

Contrary to Linearized methods, the use of a Direct-Search location is easily applicable to realistic earth models, which may show abrupt and complicated velocity variations in three-dimensions. In addition, it imposes little restriction on the form of the measure of misfit, it is stable even with poorly constrained spatial locations or origin times, and it produces comprehensive, probabilistic solutions mapping the full location uncertainty, often a complex time-space function (Lomax et al 2009 ).

This non linear probabilistic location technique has widely been used for regional hypocenter location in highly heterogeneous areas (Lomax et al. 2001 ' 200 , Zollo et al. 2002 , Husen et al. 2003 , Presti et a/. 2004 and Lippitsch et al. 2005 ), and it has been also incorporated in more local microearthquake studies (Spillmann et al. 2001 ).

We carried out all event locations with probabilistic absolute location using the Oct-Tree sampling algorithm to perform a global search within a parameter space M. This algorithm can be easily applied with 3D velocity models because it does not require partial derivate information, which is difficult or impossible to obtain in complex models (Lomax et al. 2001 ). It also permits a systematic coverage of search region and an accurate recovery of complex pdf (Lomax et al. 2009 ).

Thus, in our study the hypocenter with the maximum likelihood is determined by using the Oct-Tree nonlinear method, which consists on a recursive subdivision and sampling of rectangular cells in three-dimensional space to generate a cascade structure of sampled cells, such that the spatial density of sampled cells follows the targeted pdf values. The relative probability that a seismic source location lies within any given cell i is approximately:

(3) where V t is the cell volume and is the vector of

)

coordinates of the cell center (Lomax et al. 2009 ).

The likelihood function in Equation 1 is entirely defined by the probabilistic error processes involved. Sometimes it is desirable to change the approximations employed in Equation 2 in order to remove biases or instability in solution. Oct-tree sampling can use the L2norm misfit function or the EDT (equal differential-time) formulation (Lomax 2005 ) to represent the pdf of the data error variation, since both require searching over threedimensional spatial locations only.

We noted that for our studies the L2-norm seems most appropriate to the area dimension respect to the EDT-norm, which fits well natural earthquake location and large-scale seismic networks. The L2-norm is preferred and in this case the location likelihood function is:

pdf LSL2 =kexv\- ZJ (4)
with x -the position of the source, T obs and T cal -the arrival times observed and calculated from observation, 0 obSi and 0 cd are the polarisation angles observed and calculated, <r i and CTj are the uncertainties, and k is a normalisation constant.

The quality of the solution is assessed according to Lomax's recommendation (Lomax et al. 2009 ). We consider the maximum likelihood hypocenter, defined as the point in space of the maximum value of the location pdf (Lomax and Snieder 1995 ) and the corresponding origin time. The maximum likelihood (or minimum misfit) point of the complete, non-linear location pdf, is selected as an "optimal" hypocentre. The significance and uncertainty of this maximum likelihood hypocentre cannot be assessed independently of the complete solution pdf.

A relative measure of the volume of the high likelihood region of the location pdf, V pdf , is given by:

• pdf (x) y Pdf ~Pdf" -dV , ( 5
)
where pdf m ™ is the maximum value of the location pdf in M. The V df is a measure of the "extent", if the V df value

' pdf Pdf
is big the pdf tops up the volume. If the pdf is constant everywhere, the V pdf is equal to the volume of the model space.

The advantage of this measure is the handiness to implement it in the Oct-tree algorithm.

We also examined statistics related to the solutions using the half-lengths of three principal axes of a 68% confidence error ellipsoid approximation to the location pdf. It is an ellipsoidal, "Gaussian" or normal statistic approximation to the pdf, truncated at the 68% confidence level. This confidence ellipsoid is generated from the covariance matrix C of the pdf scatter sample (Press et al. 1992 ).

If the pdf was perfectly ellipsoidal, then there would be 68% probability that the hypocenter is inside the ellipsoid and the ellipsoid is centred on the maximum likelihood hypocenter. In fact, the maximum likelihood hypocenter can be outside of the ellipsoid -this is common with complicated or irregular (non-ellipsoidal) pdf shapes ). Anyway, it is well known that the Gaussian estimators and resulting confidence ellipsoid can be good indicators of the uncertainties in the location only in the case where the complete non-linear pdf has a single maximum and has an ellipsoidal form.

In addition, the location uncertainty can be nonellipsoidal (non-Gaussian) because the forward calculation involves a non-linear relationship between hypocentre location and travel times (Lomax et al. 2001 ).

Finally, another essential aspect for such location algorithms is the theoretical seismic wave travel-time dataset through a predefined velocity model enabling travel path calculations between any considered seismic source and receiver locations to determine the travel-time field, ray and paths take-off angles. There are three basic classes of methods to calculate travel-times and rays (full-waveform methods, ray methods, and eikonal and shortest-path, graph-based methods), among these the Fast-Marching Method (FMM) is chosen and used.

FMM solves the wavefront propagation problem through numerical solution of the eikonal equation for ray propagation, with repeated application of Huygen's principle while taking into account causality, that is to say 1999 that information only flows forward in time (Sethian ). This condition makes the method unconditionally stable in the presence of shadow zones, diffractions, etc., and also makes it applicable to 2D and 3D regular and irregular grids, such as complex 3D velocity model, including significant mining voids or paste backfill zones.

In addition, the FMM also permits through postprocessing of the travel-time field using finite-differences to follow the local time gradient back from the receiver to the source. This is of great importance for improving microseismic locations when using only a small number of sensors.

To make the location program efficient for complicated 3-D models, the travel times between each station and all nodes of an x, y , and z spatial grid are calculated once and then stored on disk as travel time grid files. The forward calculation during location reduces to retrieving the travel times from the grid files and forming the misfit function.

NUMERICAL TESTS

Synthetic Data

The localisation module implemented in SYTMISauto software used in this study, is based on the combination of the microseismic wave arrival times as well as the polarisation angles, in order to determine the hypocenter with the maximum likelihood.

Various numerical simulations were run to examine the robustness of the source location algorithm described above. Here, one of those numerical simulations is presented and it is carried out to assess the presence of an important mining void underlying a rather contrasted geological overburden made of different horizontal rock strata with different geophysical properties. This illustrative numerical test is drawn from the field situation of Cerville-Buissoncourt cavern collapse experiment during which induced seismicity of thousand events was recorded by a high resolution near filed array (Contrucci et al. 2011 ). In this case, both the modelled brine cavern and the overlying lithology have dimensions of the same order of amplitude to the distance and extension of the nearby microseismic array overlying the salt mine, approximately 120 meters for the cavern diameter, for the surface and subsurface microseismic array.

Two virtual hypocenters are positioned in the immediate vicinity of the 3D cavern, blocked and then localised (Figure 1) to compute theoretical arrival times, ray paths and take-off angles on 3D receivers. Take-off angles are expected to constrain efficiently any direct search event location method since the seismic monitoring array offer a good coverage in X and Y but not in Z (depth). The first event is positioned just above the cavity, in the cap rock, the second besides the cavity, in the salt rock.

Table I The misfits (mis.) indicate the differences between theoretical and calculated coordinates of the two hypothetical sources with known positions

These two positions are representative of the two main mechanisms to be detected within the real seismic swarm to be processed, i.e. upward progression of the cavity along with the failure of the low strength cap rock and lateral extension of the cavity in the salt. The input data, i.e. computed travel times and polarization angles, are then disturbed assuming data errors associated to sources of inadequacies and uncertainties. We perturbed the polarization angles introducing gradually an error from 5 to 15° and the P-wave pickings introducing an error range from 0.001 to 0.10 s.

The two events are then localised twice, using a 3D velocity model first with the presence of the 3D cavern and second without it. The known geometry of the cavern was approximated by a 150 m x 50 m cylinder (Figure 1).

The misfits (x, y and z) between initial theoretical coordinates and the back calculated ones are showed in the Table I.

The results show clearly that the x and y coordinates are significantly influenced by the presence or not of a large void in the velocity model, and also the calculated depth changes drastically for the source prepositioned on the sidewall of the cavity. The cause is the very different ray paths computed in the case a 3D velocity model with or without cavity. An example of ray path computed with taking into account the mining void and the clustering of the pdf points are showed in the Figure 2 In both cases the clusters of dots (Figure 2) are located around the points of the location and they are scattered, this is probably caused by the uncertainty in the spatial location due to the travel time calculation errors, velocity model, network geometry etc.

The quality of the determination of the hypocenters can be quantified and analysed using the high likelihood region of the location pdf and the half-lengths of the principal axes of a 68% confidence error ellipsoid fitting coarsely this volume (Table II). The pertinence of this last measure may be cautiously assessed through graphical representation, in fact to represent most completely the results of probabilistic, direct, global-search methodologies we used above all the geometrical properties of the location pdf . Generally, the near-ellipsoidal form of the location pdf shows a well constrained location, and it is true for this numerical simulation for both locations (see Figure 3).

In the present case, the main difference is for the first axis (1) of the ellipsoid, which fits the z axis. As expected, the most sensitive component to the strong perturbation introduced through the cavern full of brine is the vertical coordinate.

GEOLOGICAL SETTING AND MINING PROCESS

For the next step validation of the 3D source location algorithm, based a true dataset, we chose a mine induced microseismic swarm issued from the same site described here before. • A 9 meter thick, very stiff, competent layer of dolomite and anhydrites, called the "Beaumont Dolomite" strata, characterized by a 2.9 KN/m 3 density, a UCS close to 200 MPa and a velocity V P3 estimated to 4660 m/s.

• A 55 meters thick strata of marls intertwined with indurated 4620 m/s).

anhydrites ( V P4 = 3970 m/s, V P .

• The salt deposit ( V P6 estimated to 4030 m/s), lying between 185 to 340 meters deep.

From the feedback experience, the Beaumont dolomite layer plays a strong role in the overall stability of the overburden during the progressive enlargement of the brine 901 9

field (Laouafa et al. ). It qualitatively fulfils the function of a stiff elastic brittle frame able to support high stresses with low deformations before brutal failure (Nothnagel 2003 ). Nowadays, this collapse is looked for by the miners in order to prevent any post-mining risk in the long-term future. This stratum is usually seen as a clear physical barrier separating two geological horizons, the underlying large salt deposit and its cap rock and the overlying geological burden. The seismic fracturing and goaffmg process (generally associated to the longwall mining operation) of the 55 meters of cap rock may be considered then as a first critical stage of enlargement of the solution cavern. The following stage is supposed to be the progressive uncovering of the Beaumont Dolomite strata, showing a very weak, difficult-to-detect elastic response to most monitoring surface and subsurface sensors. Once the brine cavern reaches a second critical size, the Beaumont Dolomite rock strata is expected to fail brutally, leading straight to the catastrophic collapse up to the surface. In 2007, the geometry of the brine cavern was surveyed by three dimensional sonar measurements. Seismic activity was almost inexistent. In March 2008, a microseismic crisis occurred, more than 10 000 events recorded in a few days (Contrucciefa/. 2011 ).

Although there was no doubt that this microseismicity was associated to the failure of the cap rock the question raised how much of it has failed and whether or not the first phase was eventually achieved.

Fortunately, the geometry of the brine cavern was surveyed by three dimensional sonar measurements once more. This revealed that 35 meters thick by 200 meters long of cap rock strata has failed, quite asymmetrically on the east part, with the fresh top of the cavern almost licking at the Beaumont Dolomite strata. Of course, all this information may be used as a dataset for the calibration of a 3D location procedure in a complex mining environment. This means a heterogeneous rock medium hosting a large and fast changing mining void and a critical need for accurate location of the microseismic pattern to detect the entry into the last phase of development before large scale failure.

Forty microseismic events from the early stage of this seismic crisis and featuring clear easy-to-pick arrival times and polarised signals from the three components probes were selected to test and optimize the location procedure based on a 3D velocity model. Since this subset is selected on the early stage of the crisis, those events are reasonably thought to be located in the first ten meters above a non planar roof, located at an average elevation of a few meters inside the cap rock from the 3D survey. The location of those events was first computed based on arrival times and emergent angles minimizing the misfit with the L-2 norm. As regards the Cerville-Buissoncourt experiment, the orientation of the 3-D probes installed in boreholes was measured with a retrievable attitude device before grouting.

Simulations were run to test the differences between the two locations (with or without cavity into the 3D velocity model); a first localisation was done without taking into account the cavern. The results were not completely satisfying (Figure 5), because the majority of the sources are located below the cavity in the salt formation, and this does not agree with the knowledge of roof cavern evolution. The gamma ray well-logging show that massive marls roof falls of-500 000 m 3 , i.e. 4.6 Mt of material, occurred between February and May 2008 with a transverse 9008 asymmetric roof elevation (Klein et al. ), thus the majority of the sources would be located in the marl with more or less indurated anhydrites. These results lead to improve the input data for source location and above all the velocity model. Thus, the next location of these 40 events was carried out with taking into account the cavern in the 3D velocity model (Figure 6). The results are quite different from the previous. This 3D location is characterized by few events located into the cavern and the others get on towards the dolomite layer and lie into the marl with anhydrite layer. There is also a strange concentration of sources around the M 6.3 geophone (last sensor above the cylinder), which was not used to localize the sources.

Nine events are located into the clay-sandstone-marl. No events of 2008 crisis are located into the dolomite layer. In addition, three events have an aberrant location, below the cavity.

Anyway, it could be plausible that there are microseismic sources into the cavern due to insoluble blocks which fall into the brine and seismic sources into the marls with anhydrite located inside the brine, due to the roof fall.

The study area is characterized by a kind of evolution of the microseismicity that migrates in space and time and this clearly indicates a generally vertical development of the cavern dome along with an apparent lateral dissymmetry.

Concerning statistics related to the two location solutions, there is no compact location pdf in both source locations, as show in the figures 5 and 6; this irregular pdf shape is probably due to the velocity model inaccuracy (Lomax et al. 2009 ).

Another possible reason of the uncertainty source can be the absence of discontinuities in our velocity model, we introduced the different velocity values using a velocity gradient there is not a physical surface which distinguishes two different geological layers.

The ratio between the maximum axis and minimum axis of a 68% confidence error ellipsoid is always major than 1.5 (Figure 7), this determines an elongated error ellipsoid. This kind of ellipsoid highlights that the depth variation probability density is more relevant than lateral variation one and indicates above all no depth constraint (Lomax et al 2009 ). Reducing the vertical extent of the pdf requires stations at distances of the order of the source depth or less. The addition of one or more good quality S readings, especially at the closest stations, would further improve the depth constraint (Lomax et al. 2009 ).

The ellipsoid volume values are quite comparable (Figure 8). The mean difference between ellipsoid volume of 3D velocity model with cavern and without cavern is 133.54 m 3 . However, in both locations (with or without mining void in the velocity model) the weighted, root-mean square of the arrival residual (observed-calculated) times, rms, is quite small and reasonable. So, this parameter without other information is not reliable to estimate the quality of a source location.

Clearly, all these results show that there are some aspects to improve, above all related to methodology and the specific Cerville site setting. The first is to introduce into the 3D velocity model the presence of interfaces among the different geological layers. Until now the different velocity values are introduced as a kind of velocity gradient without any physical separation, i.e. between the first layer and the second layer, etc. This lack probably influences significantly the travel time paths and consequently the final results. On the contrary, for traditional ID velocity model we used defined interfaces to identify each geological layer.

Concerning the location methodology it will be also necessary to modify the assessment of statistics related to the solutions, the use of the half-lengths of three principal axes of a 68% confidence error ellipsoid is not sufficient; it indicates largely uncertain location of the hypocenter. It could be useful to adopt a classical confidence error ellipsoid with a 90% confidence error.

Finally, regarding Cerville site setting it is necessary to find a number of sources with a good quality of polarisation angles (the database with only about 40 events with polarisation angles is very small to give for the moment an exhaustive explanation of the source migration). Then, it could be interesting to improve again the velocity model, which is not a static model but evolves every time contemporary to the natural roof evolution by marl degradation in contact with the saturated brine. The 3D velocity model used derives from numerical simulation, which probably needs some new deep analysis of input parameters.

Moreover, the parameterization of the 3D V p and V s model is not based on average station spacing, the grid spacing of 5 m x 5 m x 3 m is adopted and this probably fits well with the microseismic source resolution above all in depth, but it is not adequate for the network geometry.

Nevertheless, the Cerville site is only a chance to test our algorithm.

DISCUSSIONS AND CONCLUSIONS

Nowadays, it is generally accepted that in mining environment the traditional monitoring methods cannot meet the needs of situation. In fact, the microseismic local network is considered one of the best tools to monitor.

In parallel, one effective way to improve location accuracy is to install more high resolution sensors, that is, to increase the sensitivity of the system so that more data are included in the location. Since the seismic system becomes more sensitive, it records many more events than before, which places a larger load on seismic processing staff. So, automatic source location become a fundamental aspect for mining management, but before passing to this step it is necessary to improve the source location algorithm.

Our aim is to obtain an absolute location, for microseismicity early-warning in mines, performed rapidly and in an evolutionary manner starting with the first available phase arrivals and polarization angles. An automatic and accurate source location algorithm, which gives robust and useful location information, leaves more time for seismic staff for interpretation, hazard assessment, and risk management.

The first tests to apply our source location algorithm were done using numerical simulation and then with some real data coming from the Cerville salt cavern. Based on these results, the algorithm seems to work, but as expected, our first real database is not a perfect database to study the general microseismicity of the area (it is only a part of the main crisis of spring 2008). Moreover, in this instability mechanism the inertial effect may be significant as the roughness evolution during sliding, brine infiltration and 901 9 dynamical crack propagation (Laouafa et al. ), thus the processing of small magnitude events recorded locally calls for some precautions due to some near-field conditions adverse to accuracy.

However, the gain of the 3D location (much better separation of the fracturing processes in the cap rock and on the sidewalls of the cavern, discerning unambiguously the impact of the roof falls at the bottom of the cavity) starts to be achieved and other improvements in the 3D velocity model could lead to new satisfying goals.

In addition, it is necessary to test this source location algorithm with another case of study, for instance with a rich database coming from a deep mine, this could help us to validate in a complete and evolutionary way our 3D algorithm.

Anyway, this source location algorithm is currently being improved and automated for operational implementation in a deep mine.

Figure 1 :

 1 Figure 1: Section showing the two dots (1, 2) as the hypothetical sources with known positions. The cavern is modeled by a 150 m x 50 m cylinder and the big grey spheres illustrate the seismic array

Figure 2 .

 2 Figure 2. The ray paths and the probability density of the event location are shown for the two sources re-localised with taking into account the cavern in the 3D velocity model

  The geological setting is made of different subhorizontal rock strata part of the eastern edge of the Parisian sedimentary basin, with:

Figure 3 .

 3 Figure 3. The error ellipsoid for the two sources relocalised with taking into account the cavern in the 3D velocity model

Figure 4 .

 4 Figure 4. Plan view of the instrumented site of Cerville-Buissoncourt showing microseismic stations M1-M6, the exploitation wells of line 2 100 and line 2 200, and the actual cavern limits (dashed Une, from Mercerat et ai 2010 )

Figure 5 .Figure 7 .Figure 8 .

 578 Figure5. The XY view of microseismic event location (3D velocity model without taking into account the mining void). Below an example of ray path and error ellipsoid (the small grey dots represent the probability density of the event location)

  .

		Mis.	Mis.	Mis	Mis.	Mis.	Mis.
		X	Y	Z	X	Y	Z
		(m)	(m)	(m)	(m)	(m)	(m)
		1 . At the top,	2.	On the side,
		in	the cap rock	in	the salt rock
	Cavern	0	-2.45	-7.80	5.00	7.55	19.51
	No cavern	0	-2.45	-7.80	-17.5	-14.95	14.53

Table II Example of the half-lengths of three principal axes of error ellipsoid. Axes 1,2 and 3 are semi-major (longest), the semi-intermediate and semi-minor (shortest) axes, respectively, of the confidence ellipsoid

 II

		Error Ellipsoid
		Axis 1 (m) Axis 2 (m)	Axis 3 (m)
	Cavern	69.35	37.28	36.98
	No cavern	93.52	43.84	35.02

ACKNOWLEDGMENTS

The authors would like to thank:

• SOLVAY, S.A. for the experimental site, the technical support and the data for site characterization.

• The European Commission for their support and the partners of the EU FP7 project I2Mine (http://www.i2mine.eu).

• RaSiM Organizing Committee and the reviewers for their comments that help to improve the manuscript. Special thanks also go to Anthony Lomax for his advice, counsel and expertise about the localisation methodology.


Anthony Ung
email: anthony.ung@ineris.fr

Laure Malherbe

Frédérik Meleux

Bertrand Bessagnet

Frederik Meleux

Laurence Rouil

Ensemble operational

Keywords: Air quality, modeling system, ensemble approach, validation with ground based measurement

come

Introduction

MACC-II -Monitoring Atmospheric Composition and Climate (http://gmes-atmosphere.eu)is the current COPERNICUS (previously known as GMES Global Monitoring for Environment and Security) -Atmosphere Service throughout Europe. Two operational services based on the strength of multi-model approaches are dedicated to air quality monitoring in Europe: the first one focuses on routine and near real time products (forecasts, Near Real Time analyses) and the second one on detailed analysis of the past period. Seven regional modeling teams combine state-of-the-art atmospheric models with Earth observation data assimilation. Table 1 summarizes the fixed specifications and shared input data by all models. All models use the same anthropogenic emissions dataset, the same meteorological forcing and the same lateral boundary conditions provided by the global model forecasts. This choice has been made in order to limit the spread in the model results and allow some interpretation of the differences. Nevertheless, the variability among the models is still quite large and highlight the most sensitive geographical areas where the re-analyses are still uncertain, either because of model weaknesses (uncertainties on emission, complex topography, limit of the dynamic and chemical parameterizations...) or because of a lack of available representative observations. Each team is free with the choice of data assimilation strategy (table 2) provided that the results are relevant, reliable and improve significantly the raw simulation results. The list and the current configuration of the modeling systems are detailed in QA/QC dossiers and available on the MACC project website for each model. All models have also very significant and specific skill and functionalities which allow to derived "ensemble" model estimations. This ensemble is currently based upon a median approach to obtain an average with improved skills compared to the individual models at no additional cost. Five years are available from 2007 to 2011. Figure 1 illustrates the ensemble yearly mean of the regulatory air pollutant concentrations (O 3 , PM 10 and NO 2 ) at surface level for year 2010.

Evaluation of the models and the ensemble

The evaluation and extensive analysis of the performances of the modeling systems and of the ensemble to predict the concentrations are fully documented on a yearly basis in the evaluation and assessment report [START_REF] Rouil | EVA assessment reports[END_REF](Rouil L., & 2010)). It is a systematic evaluation of the model, of the data assimilated results and of the ensemble model against a relevant set of dedicated observation data (from the European Airbase database).

MACC/EVA re-analyses are appropriate to report and comment on background air pollution (both rural and urban). The spatial resolution adopted in the current system does not allow catching local pollution situations, especially near busy roads or on industrial sites. Only observation data with "background" typology are considered.

This observation set is split into two sets: a set of stations for assimilation procedure and a set of stations for validation procedure (figure 2). Stations are randomly distributed and give a good spatial cover over European domain. The graphics illustrate the superiority of the median ensemble model (black dots) in comparison to individual models. Each model contributes to build up the skill of the ensemble median. For NO 2 , there is a limited number of available re-analysis. It results that the best model is not the ensemble but one of the member of the ensemble. Same comparison is conducted on a yearly basis since 2007 and similar results are obtained every year.

As part of the MACC-II work plan, various alternatives to the median ensemble model have been tested [START_REF] Moinat | Report on studies concerning alternative ensemble approaches, report D_R-ENS_3.2 of MACC-COPERNICUS Project[END_REF]. A wide screening of methods based at exploiting the historical knowledge of which model(s) performed "best" overall over the few last days did not give convincing results. The studies have also shown that the ensemble model takes actually benefits from the different models (including the ones with less skill), and that the "best" model, in terms of proximity to observations, is in practice changing very significantly in space and time.

From one year to another, the number and the selection of the models to compute the median can evolve, but limited to 4 or 5. As a consequence, some indicators such as the yearly mean concentration of the ensemble are dependent variables.

-RIUa is the best model for NO2 daily mean -ENSa is the best model for PM10 daily mean -ENSa is the best model for O3 daily max

Conclusion

The

Figure 1 :

 1 Figure 1: Ensemble 2010 yearly concentration for O 3 , PM 10 and NO 2 at surface level.

Figure 2 :

 2 Figure 2: assimilation (red) and validation (green) stations for NO 2 , O 3 and PM 10 pollutants. Each model performance is evaluated on the basis of classical statistical indicators which measure objectively the gap between the model results (raw data or re-analyses) and the

Figure 3 :

 3 Figure 3: Taylor diagrams for all models for PM 10 , O 3 and NO 2 with and without assimilation.

Table 1 :

 1 Fixed specification and shared input data by all models

Table 2 :

 2 Specification of the data assimilation chains run in MACC

	Characteristics		MACC2 specifications
	Domain			GEMS (Europe) Extended domain to include Turkey and Island
	Spatial Resolution		10 km
	Vertical domain (height, nb of levels)	Free
	Emission data		Improved MACC2 emission inventory + forests fire emissions
	Boundary conditions		MACC2/GRG + MACC2/AER to include global aerosol products
	Meteorological Re-analysis		IFS (ECMWF)
	Chemical Pollutants Computed		O3, NO2, PM10, PM2.5
	Chemical schemes		Free
	Data assimilation approach		Free
	in-situ data		AIRBASE (mandatory)
	Earth observations		Experiments with satellite observations
	Code	Model	Data assimilation chain	Observation data assimilated or tested in MACC
	CHM	CHIMERE	Kriging and Optimal Interpolation at ground level	Only in-situ data from AIRBASE
	EMP	EMEP	No operational data assimilation
					AIRBASE in situ measurements
					MOSAIC air borne in situ measurements
	RIU	EURAD-IM	4DVAR		NO2 tropospheric column retrievals from OMI,
					GOME-2, SCIAMACHY, MOPITT CO profiles
					PM10 SYNAER retrievals
	KNM	LOTOS-EUROS	Ensemble Kalman filter	AIRBASE in-situ data for O3 and PM10 OMI for NO2, evaluation studies
	SMH	MATCH	3Dvar with transform into spectral space	AIRBASE in-situ data Satellite information : retrieval ongoing from OMI, AIRS, MODIS, but not used operationall
					AIRBASE in-situ data
	MFM	MOCAGE	3Dvar operational chain	Satellite NO2 (OMI) and CO (IASI from ULB/LATMOS or EUMETSAT ; MOPITT from
					NCAR)
	FMI	SILAM	3Dvar operational chain	AIRBASE in situ data
	ENS	Ensemble	Median

  availability of a set of several model results allows a qualitative estimation of the model uncertainty which is linked to the variability of the model results. Southern Europe and Eastern Europe are the European sub-regions where the model responses are the most variable and uncertain, even with the data assimilation systems. This is not a surprising result; the difficulties for simulating the dynamic and chemical atmospheric processes that occur in these geographical areas are well-known from the scientific communities and are still the subject of research project. Moreover lack of measurements in Eastern Europe (compared to other geographical areas, see in annex) can also result in poor performance when modeled indicators are compared to observations. But this must be interpreted with caution.MACC-II products are considered as the "best available or most realistic" representations of air pollution patterns over Europe: more relevant than interpolation of observations and more accurate than raw simulations or one unique modeling. In its operational stage (foreseen in 2014) the GMES/COPERNICUS atmospheric services will propose complementary and comprehensive information established on the basis of state of art chemistry-transport models run operationally by modeling teams with a long experience in the field of air pollution, and extended in-situ and satellite observation sets gathered and made available. Such information could be relevant for assessing rural and urban background concentrations and exposed areas to be reported according to the air quality Directive of the European Commission (2008/50/EC).The work is done with contribution from the MACC regional modeling teams. CERFACS, FRANCE: Sébastien Massart CNRS/LISA, FRANCE: Matthias Beekman, Gilles Foret FMI, FINLAND: Mikhail Sofiev, Julius Vira KNMI, NETHERLANDS: Henk Eskes Meteo France, FRANCE : Vincent-Henri Peuch, Virginie Marécal Met.no, NORWAY: Alvaro Aldebenito, Michael Gauss FRIUKK, GERMANY: Hendrik Elbern, Elmar Friese, Achim Strunk SMHI, SWEDEN: Lennart Robertson TNO, NETHERLANDS: Arjo Segers, Lyana Curier INERIS, FRANCE : Laure Malherbe, Frederik Meleux, Bertrand Bessagnet Analyses and 3-day forecasts of European air quality are provided through the MACC-II website in near-real-time: http://www.gmes-atmosphere.eu/services/raq/raq_nrt/


Emmanuel Leprette

Stéphane COUSIN b Jef Snoeys

Christophe Proust

Improving explosion protection methods for industrial processes : the collaborative project DELFINE (Evolution des méthodes de protection des procédés industriels contre les effets des explosions : le projet DELFINE)

Keywords: Explosion, protection systems, venting, dust, turbulence

Classical methods for vent sizing are based on experimental correlations coming from academic test situations, sometime quite different from real industrial conditions. Thus, they do not take into account the variation of turbulence in industrial enclosures, which may alter drastically the explosion violence. The DELFINE installation presently in construction will allow studying dust explosions in real working conditions of a dust collector. Preliminary experiments in an instrumented small size filter are presented, which allow characterizing the turbulence level in every work phase of the system. It appears that the turbulence level is generally lower than in the standard conditions for vent testing, even during reverse jet cleaning. It should offer perspective for optimization of the venting areas for dust collectors. Explosion tests in real conditions on the DELFINE installation will be performed to try to confirm these first observations.

Background

Dust explosions continue to represent major risks in the process industries. Some accidents, such as the explosion of the grain silo in Blaye, France (August 1997, 12 victims: Masson, 1998), remind us that explosions in the industry can have serious consequences. There is an obvious need to be capable of "engineering" the safety to reduce the consequences of such accidents. In Europe, a legal framework has been imposed to implement prevention and protection measures (EU directives 94/9/CE and 99/92/CE) and practical guidelines and standards have been issued. One of these standards, EN14491 describes a comprehensive method to calculate vent areas and design venting systems. This standard is an interpretation of preceding guidelines (VDI 3673, NFPA68) by a panel of European experts. In the USA, the venting guideline NPA68 was substantially revised and the 2007 version became a standard. Vent sizing methods have received considerable attention, especially since the work of [START_REF] Donat K | Auswal und Bemessung von Druckentlastungseinrichtungen f¨r Staubexplosionen[END_REF], followed by experimental test programmes by Pineau (1982), [START_REF] Radandt S | Berufsgenossenschaft Nahrungsmittel und Gaststätten RUNES[END_REF], [START_REF] Lunn | The Effect of Vent Ducts on the reduced Explosion Pressures of vented Dust Explosions[END_REF][START_REF] Bartknecht | Explosion-schuzt : Grundlagen und Anwendung[END_REF][START_REF] Eckhoff R | Maize starch explosions in a 236 m3 experimental silo with vents in the wall[END_REF]. A number of correlations have been proposed, deriving more or less explicitly, from the theory of flames [START_REF] Rust | Explosion venting of low pressure equipment. A theoretical solution to explosion relief[END_REF]Runes, 1972). It can be recalled [START_REF] Proust Ch | Differentiated routes for the simulation of dust explosions[END_REF] that, if A is the size of the vent aperture in a vessel of volume V, the internal overpressure ΔPred ("reduced explosion pressure") can be expressed as [1] With C d the discharge coefficient of the orifice, A f the flame area, S t the burning velocity and E exp the expansion ratio of the burnt products. The "expansion velocity" S t . (E exp -1) mainly depends on the 1 emmanuel.leprette@ineris.fr ( )

2 exp 1 2 f t red d A S E P C A ρ   ⋅ ⋅ -   ∆ = ⋅   ⋅  
properties of the cloud itself (particle size and distribution, dust concentration, turbulence level, etc.) so that for a "given cloud" the explosion overpressure will correlate to the geometrical parameters implicitely contained in A f (and explicitely in A). The practical use of [1] is not straightforward because the determination of the expansion velocity S t . (E exp -1) is rather difficult. Fortunately, flame theory suggests a direct link between the expansion velocity and the maximum rate of pressure rise of the explosion in a closed vessel:

[2]

Thus, the "flame speed" parameter K ex is much more conveniently determined indirectly at standard conditions in, for instance, the 1m 3 ISO vessel (ISO6184/1) and is better known as "K st ". Note that there is no evidence that the way in which the flame propagates in these standard test conditions corresponds to any practical situation. However since the experimental conditions are kept constant in this standard testing [START_REF] Proust Ch | Measuring the violence of dust explosions with the 20 liter sphere and the standard ISO 1 m3-systematic comparison and interpretation of the discrepancies[END_REF], the variation of the K st between dusts certainly reflects the differences in terms of powder reactivity.

Usually the flame theory cannot be applied directly and some fitting with the experimental data is required. Some of those fitted correlations have been incorporated into guidelines and standards (NFPA 68, VDI 3673, EN 14491). Most process parameters suggested by [1] are covered such as the volume of the vessel, the shape factor from the proportionality between A ch (internal area of vented vessel) and V 2/3 and the reactivity of the dust cloud via K St and P max . The European standard (EN 14491) for example proposes the following relationship [3]:

where L/D represent the shape factor of the vessel. Eckhoff discussed this kind of approach and recalled that K st is not a material constant [START_REF] Eckhoff R | Dust explosion experiments in a 500 m3 silo cell[END_REF][START_REF] Eckhoff R | Maize starch explosions in a 236 m3 experimental silo with vents in the wall[END_REF][START_REF] Eckhoff R | Maize starch explosions in a 236 m3 experimental silo with vents in the wall[END_REF] and should depend also on the turbulence. Eckhoff produced data in very large silo explosions which showed that the severity of the explosion was greatly affected by the state of the dust cloud at ignition. Further controlled experiments (table 1 from the data of [START_REF] Eckhoff R | Dust explosion experiments in a 500 m3 silo cell[END_REF][START_REF] Eckhoff R | Maize starch explosions in a 236 m3 experimental silo with vents in the wall[END_REF][START_REF] Eckhoff R | Maize starch explosions in a 236 m3 experimental silo with vents in the wall[END_REF][START_REF] Tonkin | Dust explosions in a large scale cyclone plant[END_REF][START_REF] Pineau | Gas and Dust Explosions in closed and vented Vessels[END_REF], including large tests which more closely resembled real processes, facilitated the confrontation between the "K st " based correlation (as [3]) and the realistic experimental results.

A comparison between calculated vent areas (by using [3]) and experiments in close-to-real configurations is shown in figure 1. The calculated vent areas diverge significantly from the measurements with some significant over estimations of ∆P red but also severe underestimations.

Figure 1: Calculated (by using [3]) and measured overpressures for vented dust explosions developing in real process equipment

This situation is by no means acceptable from a safety point of view and not even more from the industrial point of view because the safety performance is not ensured. One likely reason for these discrepancies is the significant influence of the initial turbulence of the cloud [START_REF] Amyotte | A review of the effects of turbulence on gas phase, dust and hybrid mixture explosions[END_REF][START_REF] Eckhoff R | Dust explosion experiments in a 500 m3 silo cell[END_REF][START_REF] Zalosh R | de petite taille sont présentées, qui permettent de caractériser le développement de la turbulence dans les différentes phases de fonctionnement du système[END_REF][START_REF] Proust Ch | Differentiated routes for the simulation of dust explosions[END_REF] on the combustion (via S t ). Because the turbulence level depends strongly on the flow field inside the vessel prior to ignition, it is natural to think that a more accurate prediction of the course of explosions can only be achieved if sufficient details on the flow field could be entered into the combustion model (incorporating explicitly the influence of turbulence). One of the aims of DELFINE project is to be capable of better predicting the course of the explosion to achieve a better explosion projection by venting.

DELFINE experimental setup

Experiments will firstly be performed in filtration process in a filter not only because dust collectors are involved in a majority of the reported dust explosions, but also because the reverse jet cleaning system (RJC) -very frequently used in dust collectors-is often invoked as being the reason for explosions.

The experimental setup will consist in a DELTA NEU dust collector (volume about 4 m 3 ) connected to an open dust circulation duct (figure 2). The duct will be equipped with a dust injection system (constant mass flow rate injector). The injector is designed to deliver dust fuel to a flowing system in a known feed rate, thereby allowing for the creation of a homogeneous concentration atmosphere in the whole system.

The dust collector will be modified for the implementation of various explosion protection systems (venting, suppression, isolation) at different locations, and will offer many ports with large possibilities for measurement: turbulence, concentration, flame and pressure. Additionally videos from the inside of the collector and the pipe will be captured through transparent walls and windows.

Figure 2: View of the future DELFINE experimental set-up

As a first step, preliminary work was performed in a real dust collector at small scale to test the instrumentation and help the design of the future installation. A 1/3 scaled down filter + pipe device available in the premises of one the partners (FIKE) was used to prepare the instrumentation and begin the physical analysis. A full operational dust collector with piping was purposely built (figures 3 and 4). Two different lengths of filter bags or 'envelopes' are used to represent a small (index 9 on figure 4) and larger filter area (index 8), and different distances from filter bag to filter bag and filter bag to walls. The flow rate can be adjusted. The reverse jet cleaning system consists of a 1 L tank pressurized to 5 bar, solenoid valves release short (typical duration 0,2 sec) blasts of air into the filterbags (index A). During RJC, short blasts of high pressure air are introduced through venturi's and nozzles inside the sleeves (figure 5). This pulse of air travels through the filter sleeve, causing it to flex out from the cage to dislodge the dust accumulated on the filter sleeve (figure 6). During this period of time, the flow is globally "reversed", large amounts of dust fall into the hopper and it is believed that a significant degree of turbulence appears. In the real system the turbulence should be measured while a large quantity of dust will be conveyed so that LDA techniques will not work anymore. To overcome this very important difficulty, the bidirectional probe system (McCaffrey, 1976) once used by Tamanini (1990) to measure turbulence was further developed (figure 7). The gauge is a small portion of a cylinder with a separating wall in the middle (10 mm outer diameter, 20 mm long). The difference of pressure measured on both sides of the separation is proportional to the dynamic pressure applied along the axis.

Figure 7: Schematic principle of bidirectional probe

The instantaneous velocity is deduced directly from the measurement. A Furnace Control type of sensor was used enabling velocity measurement in the range ±10m/s with a cutoff frequency of some tens of Hz. having a sufficient dynamics (10 Hz) and accuracy (±100 Pa). A detailed description of the techniques falls out of the scope of this paper but the results obtained with this technique can be compared with known data. For instance, the mean velocity (U) and turbulence intensity (u' : root means square of the fluctuations of the velocity around the mean value) of the flow was measured inside the feeding pipe (at a location where the flow field was established, far from the inlet) and at several position across the section. An excellent agreement is found between the measurements and the available data from the literature (figure 8).

Flowfield and turbulence inside the experimental setup

Measurements were performed at several locations inside the filter (figure 9), during normal steady state operations and during reverse jet cleaning. All the tests were performed without dust added although it is possible to do the measurement with a dust cloud with the bidirectional probes.

Steady state operation

The flow in the inlet pipe (location 6) during steady state operation was measured as outlined above described above (figure 10 and 11): with an average velocity of 7 m/s (resp. 10 m/s), the turbulence intensity was on the order of 0.4 m/s (resp. 0.6 m/s).

Location 1 represents what happen just below the bags in the empty space of the filter. It can be seen (figure 10) that gauges 1&2 detect a globally positive velocity2 , gauges 4&5 see a negative flow whereas for the central gauge (n°3) the average velocity seems close to zero. This suggests a swirling flow. The average velocity over all the probes amounts about 2 m/s (absolute value) when the flow velocity in the pipe is set to 10 m/s. Since the cross section of the filter is about 0.16 m 2 and that of the pipe 0.03 m 2 , the average velocity seem to be in proportion to the cross section ratio. Note however that the turbulence intensity is significantly larger than in the pipe, typically more than 1.1 m/s, which might result from the velocity gradient induced by the swirling flow. In location 2 (between the sleeves), the flow seems mainly orientated in the same direction with an average velocity (calculated on all the gauges) of about 1m/s (absolute value). The turbulence intensity is about 0.8 m/s, smaller that in location 1. The presence of the sleeves may dampen slightly the velocity fluctuations. A summary of all the measurements done at all locations is presented in figure 11.

During reverse jet cleaning operation

A significant influence of RJC operation on the turbulence is expected. Statistical averaging was applied on a batch of about 30 identical pulses (figure 12 for two gauges at location 1). Note that the average flow is largely positive so that it is now going from the top toward the bottom of the filter with a significant velocity, peaking up at about 8-10 m/s during the pulse. It is on the same order than the steady state flow at this particular location. Fluctuations in time are noticeable, with a frequency of 15 to 20 Hz which seems a natural acoustic resonance of the filter/duct assembly. This last remark suggests that rather than producing a turbulent flow, the air-pulse, triggers the acoustic vibration of the filter with is more a kind of plug pulsating flow with no randomly distributed eddies (which we would expect from a turbulent flowfield). This result would suggest that this kind of RJC process would not degenerate in more turbulence than in the mean flow during normal filtering operations. This important point will be verified in the final DELFINE setup using transparent walls.

Implications and conclusion

In vent sizing equations, the severity of the explosion is represented by the coefficient K St . It is measured in very specific turbulence conditions (ISO6184-1) which were characterised in a preceding paper [START_REF] Proust Ch | Measuring the violence of dust explosions with the 20 liter sphere and the standard ISO 1 m3-systematic comparison and interpretation of the discrepancies[END_REF]. In particular, the turbulence intensity in the ISO chamber is about 2 m/s which is about twice that measured in the filter. In another paper (Figure 13 from [START_REF] Schneider | Determination of turbulent burning velocities of dust air mixtures with the open tube method[END_REF], a link could be made between u' and S t which should be roughly proportional to K st (expression [2]). It can be realized that S t , and following Kst, is roughly halved when u' drops from 2 m/s to 1 m/s. Refering to expression [3], all other parameter being conserver, a decrease of Kst by 50% would result in a decrease of 50% for the required vent area.

Figure 13: Turbulent burning velocities measured in a 30 cm diameter tube (starch-air clouds from [START_REF] Schneider | Determination of turbulent burning velocities of dust air mixtures with the open tube method[END_REF]

Figure 3 :Figure 6 :

 36 Figure 3: Purposely built dust collector

Figure 8 :

 8 Figure 8: arrangement of probes in the inlet tube and results (curves : original graphs from Comte-Bellot and Hinze stars : present measurements) no dust in the flowfield

Figure 9 :

 9 Figure 9: locations of the measuring points in the filter (location 6 is on the inlet pipe, 1&3 below the sleeves, 2 between the sleeves, 4 top between wall and short sleeves)

Figure 10 :

 10 Figure 10: flowfield at location 1 (left: inlet velocity10 m/s) and 2 (right: same velocity)

Figure 12 :

 12 Figure 12: RMS velocity and flow velocity modulus-location 1 (during reverse jet cleaning pulse)

Table 1

 1

		: vented dust explosion experiments in real or "realistic" equipments
	ref Author	Equipment Dispersion of the powder	Injection	Pipe	Volume	Aspect	Area of	Dust
					velocity-	diameter	of the	ratio	the
					pressure	(m)	vessel		opening
					(m/s-barg)		(m3)		(m2)
	A Eckhoff,	filter	tangential flow of dust from a	35 -0	0,155	5,8	2	0,11	Maize
		1988		pneumatic transport line						starch
	B	Eckhoff,	filter	tangential flow of dust from a	35 -0	0,155	5,8	2	0,2	Maize
		1988		pneumatic transport line						starch
	C	Eckhoff,	filter	tangential flow of dust from a	35 -0	0,155	5,8	2	0,4	Maize
		1988		pneumatic transport line						starch
	D Eckhoff,	filter	tangential flow of dust from a	35 -0	0,155	5,8	2	0,55	Maize
		1988		pneumatic transport line						starch
	E	Tonkin,	cyclone	tangential flow of dust from a	13 -0	0,23	1,2	1	0,06	Wheat
		1972		pneumatic transport line						flour
	F	Tonkin,	cyclone	tangential flow of dust from a	13 -0	0,23	1,2	1	0,08	Wheat
		1972		pneumatic transport line						flour
	G Tonkin,	cyclone	tangential flow of dust from a	13 -0	0,23	1,2	1	0,12	Wheat
		1972		pneumatic transport line						flour
	H Eckhoff,	silo cell	coaxial flow of dust from a	38 -0	0,155	236	6	3,4	Maize
		1986		pneumatic transport line						starch
	I	Eckhoff,	silo cell	coaxial flow of dust from a	38 -0	0,155	236	6	5,7	Maize
		1986		pneumatic transport line						starch
	J	Eckhoff,	silo cell	coaxial flow of dust from a	12 -0	0,2	500	3,5	2	Wheat
		1984		pneumatic transport line						dust
	K Eckhoff,	silo cell	coaxial flow of dust from a	12 -0	0,2	500	3,5	3	Wheat
		1984		pneumatic transport line						dust
	L	Eckhoff,	silo cell	coaxial flow of dust from a	12 -0	0,2	500	3,5	5	Wheat
		1984		pneumatic transport line						dust
	M Eckhoff,	silo cell	coaxial flow of dust from a	12 -0	0,2	500	3,5	8,8	Wheat
		1984		pneumatic transport line						dust

The sign depends on the orientation of the probe. Positive here means from the top to the bottom of the filter