Mercerat Diego
email: diego.mercerat@mines.inpl-nancy.fr

Souley Mountaka
email: mountaka.souley@ineris.fr

Bernard Pascal
email: bernard@ipgp.jussieu.fr

INDUCED SEISMICITY IN A SALT MINE ENVIRONMENT EVALUATED BY A COUPLED CONTINUUM-DISCRETE MODELLING

Keywords: Sismicité induite, cavité saline, banc raide, in situ, modélisation continue-discrete, microfissures

ABSTRACT: With the objective to better understand the induced microseismicity in a salt mine environment due to an underground solution mining, an in situ experiment is undertaken by GISOS 1 in the Lorraine salt basin. The overburden overlying the salt cavity is characterized by the presence of two competent layers where most microseismic events are expected. This paper presents a coupled continuum-discrete modelling to simulate the mechanics of fracture initiation and propagation in the rock mass overlying the cavity: discrete approach for the competent layers and continuum approach for marls, salt and other rocks and soils. For the competent layers, numerous calibrations of the model microparameters based on the laboratory results are firstly performed. The first coupled modelling results suggest that the mechanism of fracturing in the competent layers is predominantly tensile as it could be expected. The results also show that the microseismic events associated to the progressive damage in the competent layers through microcraks development can be modelled. This opens interesting perspectives to assess the feasibility of seismic monitoring of underground cavities by comparing, in the future, the numerical modelling results with the recorded seismicity of the study area. KEYWORDS: Induced seismicity, salt cavity, competent layer, in situ, continuum-discrete modelling, , microcraks, RESUME : Dans l'objectif de mieux comprendre la microsismicité induite dans des mines de sel, une expérimentation in situ a été entreprise par le GISOS sur une cavité saline de la région Lorraine. Le recouvrement est caractérisé par la présence de deux bancs raides dans lesquels l'essentiel de l'activité microsismique est attendue lors de la reprise de l'exploitation. Ce papier présente une modélisation basée sur une approche couplée : continue (pour les marnes et le sel) et discrète (pour les bancs raides) pour évaluer les mécanismes de développement des fractures dans le recouvrement. De nombreuses calibrations des microparamètres du modèle discret à partir des paramètres macroscopiques ont été nécessaires. Les résultats de la modélisation montrent que le mécanisme de fracturation dans les bancs compétents résulte essentiellement d'efforts en traction. Les premiers résultats montrent aussi que les évènements microsismiques peuvent très bien être mesurés à travers la formation des microfissures. Ceci ouvre des perspectives intéressantes pour des confrontations avec les mesures qui seront enregistrées lors de l'effondrement provoqué.

Introduction

A research project within the GISOS program has been launched in order to assess the feasibility of seismic monitoring of underground growing cavities due to salt dissolution (solution mining). Within such framework, this specific project focus on two main complementary axes : (a) the validation of induced seismicity monitoring techniques in salt mine environments, and (b) the numerical modelling of deformation and failure mechanisms with their associated acoustic emissions as well as the induced microseismicity.

Previous studies on induced seismicity of underground growing cavities [START_REF] Johnston | A survey of mining associated seismicity[END_REF]Einstein 1990, Mendecki et al. 1999) confirm the presence of microseismic events associated with the progressive damage of the overburden until its final collapse. Continuous monitoring of such microseismicity may provide crucial information for stability analysis of the area surrounding the cavity. In order to validate such seismic monitoring techniques, a microseismic network has been installed over a site test at Cerville-Buissoncourt (Lorraine, France), where a stable brine cavity is nowadays found. When the salt exploitation restarts, the spatial and temporal evolution of the seismic event locations can firstly be recorded and secondly be correlated with the results of mechanical modelling.

Numerical modelling of such complex process requires a large scale model which takes into account both the cavity evolution within the salt layer and the mechanical behaviour of the overburden layers where high deformation and fracturing are expected. To achieve this scientific challenge (experimental and numerical through geophysics and geomechanics fields), a coupled continuum and discrete modelling strategy is then adopted. This paper focus on the feasibility to achieve the previous scientific objective and presents the preliminary numerical results of induced microseismicity evolution in the competent layers through the microcracks distribution.

From numerical point of view and in order to take into account the elasto-brittle behaviour of the competent layers where most seismic damage is really expected, we use the code PFC2D (Particle Flow Code in two dimensions) based on the discrete element method. To describe the behaviour of the other layers (mainly composed of marls and salt) which present more ductile and/or viscoplastic behaviour, a continuum approach based on the bidimensional Fast Lagrangian Analysis of Continua code (FLAC) is used allowing for different and more complex constitutive laws.

Numerous calibration process are firstly performed to estimate the microparameters used in PFC2D to reproduce the macroscopic stress-strain curves obtained from laboratory tests (macroscopic behaviour from uniaxial and triaxial tests).

A large scale continuum model assuming an elastic perfectly plastic behaviour of the competent layers is firstly built in order to establish the suitable boundary conditions with regard to the brine cavity sizes. Finally, a fully coupled continuum-discrete model is developed (coupling between FLAC and PFC2D domains is achieved through a special function written in the internal language of both codes), as well as a comparison with the previous continuum approach.

Site test of Cerville-Buissoncourt

Brine cavity and overburden characteristics

The brine cavity of Cerville-Buissoncourt (Lorraine, France) is actually located within a salt layer at about 180 m depth and it presents an irregular shape but in a first approximation, it can be considered as a cylindrical cavity of 100 m of diameter and 50 m high. An schematic geological section of the study area is shown in Figure 1. Typically, the overburden is characterized by the presence of a series of intercalated anhydritic marls lying over the salt layer, and on top of them, the presence of two competent layers of 8.5 m thickness with an elasto-brittle behaviour located approximately 50 m above the salt layer. These layers correspond to what is known as the Beaumont dolomite (a dolomite of 5 m thick and an anhydrite of 3.5 m thick).

Microseismic network

A microseismic network has been installed over the underground brine filled cavity of Cerville-Buissoncourt. The network includes four 3D components and three 1D component seismometers located at different depths in cemented boreholes distributed around the study area. A velocity model of the area, from the processing of a 2D seismic line (work in progress), will allow us to obtain an accurate location of microseismic events when the salt exploitation restarts . The spatial and temporal evolution of the seismic event locations can then be correlated with mechanical modelling results, and therefore a quantitative analysis of seismic and non-seismic deformation will be available.

Geomechanical laboratory characterization

A large number of laboratory tests have been carried out on SOLVAY S.A samples to characterize the geological formations of interest. Some of them are available (by courtesy) for this study. In particular, the Beaumont dolomite and the anhydrite facies present a rather high Young modulus (~80 GPa) and uniaxial compressive strength (UCS ~200 MPa) for typical sedimentary rocks. The tensile strength measured by Brazilian tests is also high. The acoustic properties of these competent layers are quite elevated.

On the other hand, the mean values found for samples extracted from the marls located between the salt layer and the Beaumont dolomite are much lower, with extremely low tensile strengths measured in some samples.

The results of this geomechanical laboratory characterization are used to choose the physical and mechanical parameters of the numerical model, and to obtain the appropriate microproperties, by means of a calibration process that will characterize the PFC2D inclusion within the continuum. For numerical modelling purposes, we do not differentiate between the two competent layers (Beaumont dolomite and anhydrite facies).

4.

Modelling strategy

Philosophy and coupling continuum-discrete modelling

When the salt exploitation restarts, the cavity will progressively grow leading to develop irreversible damage within the upper layers until the final collapse. Numerical modelling of such a complex process requires a large scale model which takes into account both the cavity evolution within the salt layer and the mechanical behaviour of the overburden layers where high deformation and fracturing are expected.

From previous mechanical laboratory investigations, both competent layers exhibit a clear elastic and brittle behaviour. To keep the elasto-brittle behaviour of the competent layers where most of seismic damage is expected and in order to directly compare the acoustic emission (AE) events, we use the PFC2D code. That allows us to represent the damage within the competent layers directly by the formation and tracking of microcracks. A detailed description of the PFC2D model and the theory is provided elsewhere (see for example [START_REF] Hazzard | Micromechanical modelling of cracking and failure in brittle rocks[END_REF][START_REF] Potyondy | A bonded particle model for rock[END_REF]).

Below we briefly present the main characteristics of this code.

To describe the behaviour of the other layers (mainly consisted of marls and salt) which present more ductile and/or viscoplastic behaviour, a continuum approach based on the code FLAC (a 2 dimensional finite difference code widely used by the scientific and industrial communities) is employed. This leads to approach more complex rheological behaviours based on continuum mechanics (elastoplastic with or without hardening and softening, viscoelastic or viscoplastic). Hence, the brittle layers are modelled in PFC2D as an inclusion of bonded circular particles embedded in a classical FLAC grid, the applied forces received from PFC2D act on the grid points along the inclusion boundary. The translational velocities of these boundary particles are received from FLAC. The full coupling between FLAC and PFC2D boundaries is achieved through a special function written in the internal FISH language.

Short description of PFC2D code

The PFC2D code models the rock mass as an assembly of particles bonded together to reproduce cohesion. The bond that connects two neighboring particles can break apart if the normal or shear stress at the bond exceeds its normal or shear strength, respectively. This corresponds to the explicit formation of a microcrack, with the spontaneous generation of a stress wave that can propagate outward from the crack. Therefore, PFC2D was extensively used to model acoustic emission in rock specimens (Hazzard andYoung, 2002, 2002b) and to simulate fault zones (Hazzard et al., 2002). Interesting studies were conducted to associate many cracks that occur very close in space and time, into one microseismic event of spatial dimension of several broken bonds (Hazzard and Young, 2002b). This allows to explicitly assign seismic source parameters, such as moment tensor or magnitude values, emergent from the modelled material that can be compare with those obtained from observations. However, in our numerical simulations, not all zones in the numerical model are necessarily submitted to large stress variations and exhibit contact breaks. Apart from that, despite the extensive use of PFC to model elasto-brittle failure in rock and some other interesting progress (Dedecker et al., 2004), it is not clear at the present time how to model more complex rheological behaviours in PFC2D. Only the competent layers overhanging the salt cavity will be modelled with this code.

PFC2D microproperties calibration

In PFC2D, the particle contact stiffness, the contact coefficient of friction and the contact bond model between different particles govern the mechanical behaviour of the material. These are known as the microparameters (or microproperties) of the bonded particle model. They must be chosen in order to reproduce the relevant material properties measured in the laboratory tests (at a macroscale) of the rock specimen to be modelled. This calibration process is not specified a priori and a trial-and-error procedure is generally used. Nevertheless there is some knowledge and guidance based on previous research on PFC2D calibration [START_REF] Kulatilake | Physical and particle flow modelling of jointed rock block behavior under uniaxial loading[END_REF][START_REF] Potyondy | A bonded particle model for rock[END_REF]) about the relation between each microparameter and its effect on mechanical properties (Table 1).

Following this guidance and using the previous laboratory results, we have calibrated our discrete model of the Beaumont dolomite by simulating a Brazilian test and three biaxial compression tests (at 0.1, 2 and 5 MPa of confining pressure). We used a parallel bonded model material because it seems to well approximate the mechanical behaviour of a brittle elastic cement joining two bonded particles.

This calibration process is firstly performed on the more realistic particle arrangement; say a random packing arrangements around a mean grain size. However, when the specimens are created using a random packing of particles, an extra effort is needed to obtain a uniform stress state within the sample with low 'locked-in' stresses, and to eliminate 'floaters' (particles with less than three contacts with neighbouring particles). For this reason, the calibration process is also performed on a regular hexagonal packing of particles. Some differences in the microparameters and tensile-shear microcrack initiation thresholds are therefore to be expected, particularly for microcrack densities which are dependent on the porosity of the particle packing arrangement. We have firstly verified the independence of the mechanical properties to the mean particle radius by performing ten different biaxial test simulations with similar microproperties. For each mean particle radius, we only change the random number generator. As expected, when increasing the resolution (increasing the number of particles along the shortest dimension of the model), the coefficient of variation (ratio of standard deviation to mean) of each material property decreases as it was confirmed by [START_REF] Potyondy | A bonded particle model for rock[END_REF]. From these tests, the retained microparameters are given in Table 2 and the emergent macroproperties are shown in Table 3. The calibration of a regular hexagonal packing model was carried out starting from the retained values of the previous model shown in Table 2 for Davg (average particle diameter) of 0.5 cm. This regular arrangement of particles allows a straightforward installation of a required uniform stress within the sample by slightly modifying the particle radius. This may be useful for later applications regarding the coupling with FLAC code.

Particles

As expected, the stiffness and unconfined compressive strength of a regular packing specimen are much higher than those corresponding to a random packing specimen. For this reason, in order to keep an identical response, we drastically decrease both the bond strengths and stiffnesses of the regular hexagonal model (compare Tables 2 and4). Similar behaviour in terms of independence of microparameters with particle radius and upscaling to regional dimensions are obtained. In Table 4 we show the retained microparameters for the hexagonal packing bonded particle model that reproduce the mean elasto-brittle behaviour of the Beaumont dolomite.

Geomechanical model

The geometrical model shown in Figure 3 is essentially composed of :

-the salt formation extending from the bottom boundary of the model up to the actual brine cavity roof, -on top of the salt layer, a formation of intercalated marls composed of thin anhydritic to clayley marls with relatively weak mechanical characteristics, -and the overburden characterized by two competent layers (dolomite and anhydrite level), followed by a series of carbonates, sandstones and marls layers up to the surface.

The geomechanical characteristics of the salt, the competent layers and the intercalated marls are known from the laboratory tests whereas the characteristics of the other materials are inspired on the literature of these basins. For the PFC2D inclusion (in the coupled scheme) we use the microproperties previously calibrated. The inclusion is also represented in the model. In a first time and in order to calibrate the boundary conditions, many simulations have been carried out using an elasto-plastic behaviour for the competent layers, the corresponding optimal geomechanical model is then used in this study. The final geometrical model is 4600 m large and 850 m high.

Because of the absence of symmetry of the cavity, the problem is thus three-dimensional. In spite of its three-dimensional nature and for reasons of simplicity, we chose to perform our analysis in 2D for these first calculations. While considering a section in the greatest dimension of the cavity, we place ourselves in conditions of plane deformations. Consequently, the boundary conditions are:

-no vertical displacement at the bottom boundary of the model (assumption of fixed substratum).

-no horizontal displacement at the lateral boundaries (they are placed far enough to neglect any influence from the cavity filling out process).

The initial stress state is assumed to be isotropic: this assumption is well justified within the salt layer where null initial deviatoric stress is expected. In absence of data concerning the initial stresses in the other materials, we also admit that they are isotropic (a sensitive analysis on the stress anisotropy ratio in such formations can be undertaken on the basis of the literature data). Finally, the vertical stress is equal to the overburden weight.

The initial conditions within the PFC inclusion and in the cavity before the exploitation restarts:

-we established an unifrom isotropic stress of 3 MPa within the PFC domain corresponding to the approximate initial stress found at the inclusion depth due to the overburden weight .

-When the cavity is full of brine, a normal stress acting on the cavity walls equal to λ σ 0 is imposed (σ 0 is the overburden weight at the corresponding depth, and λ=0.5). The progressive filling out of the cavity is modelled by progressively reducing λ from 0.5 to 0 (empty state of the cavity).

5.

Results of large scale modelling

Fully continuum approach

As previously mentioned in the continuum approach, we assume that all off the competent layers are viewed as an elastoplastic material. The aim of this assumption is: (a) to guarantee the type of boundary conditions (applied forces or fixed normal displacements) that will be specified, as well as the external model dimensions with regard to the cavity size, and (b) to keep the model reaction (global response under the initial conditions and when filling out the cavity).

The vertical displacement evolution on top of the Beaumont dolomite and the subsidence at the topography surface are illustrated in Figure 4 In Figure 5, we can appreciate the state of the plastic zones through the layers lying on top of the cavity for λ around of 0.1. The competent layers remain stable and the plasticity migrates from the weakest layers within the intercalated marls and located at the immediate roof, up to the base of the competent layers. For λ = 0 , most rock materials intercalated between the salt formation and the dolomite fail, leading to a significant stress report in the competent beds.

Coupled continuum-discrete approach

The coupling scheme embeds a PFC2D inclusion region within the FLAC grid (the complementary domain of the inclusion in the competent beds where microseismic events are not expected is modelled as a continuum domain). A PFC2D material must be created to confirm the inclusion region such that a single layer of particles overlaps the boundary as illustrated in Figure 6. The applied forces received from PFC2D act on the grid points along the inclusion boundary. The translational velocities of these boundary particles are received from FLAC after one step of calculation (i.e. after the motion calculation where velocities and displacements are computed based on the internal and external forces acting on the gridpoints).

Using the calibrated hexagonal parallel bonded model 2 , a rectangular PFC2D inclusion within the FLAC model has been created. As a first step, we decided to limit the inclusion to the region where the greatest variation in ambient stress are produced due to the cavity growing. The inclusion centred on top of the actual cavity, is 260 m long by 8.5 m thick, and it covers the dolomite and the anhydrite levels (see Figure 3). Filling this area with circular particles of 0.5 m diameter gives a total number of 10800 particles forming the inclusion. A detail of the coupled model can be seen in Figure 6.

Before establishing the coupling process, a uniform isotropic compressive stress of 3 MPa is obtained within the PFC2D inclusion prior to its connection to the FLAC grid. This allows faster convergence of the initial coupled model to be stabilized. For the rest, the coupled FLAC-PFC2D

2 Due to difficulties encountered in installing an initial stress state near 3 MPa with a random particle packing, we have use for this feasibility study the regular particle packing, some differences are therefore to be expected.

scheme is run with similar boundary and initial conditions as the fully FLAC simulations specified in the previous section. The resulting force distribution in the cement between particles of the PFC inclusion can be seen in Figure 7. From λ = 0.5 down , we can appreciate the increasing zone in tension at the middle of the inclusion, where almost all microcraking occurs. A number of tensile cracks develop first with a random spatial distribution and concentrating later in the central part of the inclusion (just on top of the cavity). The number of microcraks (also shown in Figure 7) increases from about 20 at λ = 0.5, up to more than 100 at λ = 0.1. There is a sudden rise up between λ = 0.3 and λ = 0.1 manifesting the failure instability within the stiff layer. However, a clear tensile fracture can not be identified at this stage, as it was the case for the fully continuum modelling approach (see Figure 5).

Figure 1 .

 1 Figure 1. Schematic section of the salt cavity and main geological facies.

  k n /2 = 65 GPa σc = 120 ± 60 MPa k n /k s = 2.5 μ = 0.5 τc = 240 ± 60 MPa

Finally, Figure 2

 2 (a) compares the simulations of uniaxial test with random and regular particle packing and the observed ones. Figure2(b) illustrates the development of microcracks at different stage of axial stress (pre-peak corresponds to 50% of peak; peak and post-peak related to the end of test).

Table 4 .

 4 Model microproperties for the competent layer samples. Hexagonal packing. σc = 25 ± 10 MPa k n /k s = 2.5 μ = 0.5 τc = 50 ± 10 MPa

Figure 2 .

 2 Figure 2. Uniaxial compressive test : (a) Example of axial stress-strain curve for regular and random packing and the experimental curve, (b) Microcracks distribution (A random packing, B regular hexagonal packing) at three different stages of axial stress (tensile cracks in red, shear cracks in black)

Figure 3 -

 3 Figure 3 -Detail of the geomechanical model used in this study: the cavity in the salt layer and the PFC inclusion (used later in the coupled scheme) within the competent layers are shown.

  as a function of the brine pressure within the cavity: (a) λ = 1 virgin rock mass, (b) λ = 0.5 initial state of the cavity full of brine, (c) λ < 0.5 simulating a progressive filling out of the cavity (λ = 0 representing the empty state of cavity).

Figure 4 .

 4 Figure 4. Vertical displacement evolution at the surface and at the top of dolomite (the salt cavity is schematically shown).

Figure 5 .

 5 Figure 5. State of plasticity in the overburden layers for λ = 0.1. The competent layers remain stable (without plastic points).

Figure 6 .

 6 Figure 6. Detail of the right limit of PFC2D inclusion within the FLAC grid. Boundary particles (red) and grid node positions (black) are shown in the PFC2D domain.

Figure 7 .

 7 Figure 7. Bond forces in the cement between particles (black = compression, red = in tension) of the PFC inclusion at (from top to bottom) λ = 0.5, 0.4, 0.3, 0.2 and 0.1. Tensile cracks are shown in blue.

Table 1 .

 1 Relations between microparameters and mechanical properties.

	Microproperties	Macro mechanical property
	Normal (k n ) and shear (k s ) contact stiffness	Young's modulus (E) and Poisson's ratio (ν)
	Normal (σc) and shear (τc) bond strengths	Uniaxial Compression Strength (UCS)
	Friction coefficient (μ) -Particle clustering Ratio shear/normal bond strength	Slope of failure envelope
	Mean bond strength (σc) -Particle radius	Brazilian strength (σ t )
	4.3.1 Random packing

Table 2 .

 2 Model microproperties for the competent layers samples (parallel bonded model). Random packing.

Table 3 .

 3 Macroproperties of the PFC2D model of the competent layers samples (80 x 160 mm 2 specimens).

	Property	PFC2D model D avg = 0.5 cm
		(number of samples = 10 )
	E (GPa)	78.3 ± 2.4
	ν	0.261 ± 0.031
	UCS (MPa)	173.2 ± 8.6
	φ (deg)	42 ± 4
	c (MPa)	37.9 ± 5.1
	σ t (MPa)	44.1 ± 4.5
	4.3.2 Regular hexagonal packing

ACKNOWLEDGEMENTS : This works has been performed within the GISOS framework. We acknowledge with gratitude the financial support of the MEDD and MinEFI French ministries. to subventions. The authors also thank the company SOLVAY S.A. for the site, the technical support at the time of the experiments and the data of the site characterization

6.

Concluding remarks

We have presented a preliminary investigation of a coupled continuum-discrete modelling strategy to analyse a complex geomechanical and geophysical process in a salt mine-environment. This involves the fracturing and deformation of the overburden due to brine pumping out and the consequent pressure decrease in the underground cavity.

Based on the results of a laboratory geomechanical characterization, we have calibrated the microproperties of the PFC2D model to be used as the discrete inclusion within the competent layers in the overburden. We used both a random and hexagonal particle packings to create the parallel bonded specimens. The influence of the particle size and the scale effect was analysed. We do not notice any drastic change in the emergent mechanical properties of the modelled materials at both the laboratory and regional scale. This suggests that the use of PFC for modelling the elasobrittle behaviour of competent layers at the regional scale should be satisfactory.

On the other hand, the effect of particle packing seems to be extremely important in the macro response of the modelled material. The samples created with an hexagonal regular packing are stiffer than those created with a random packing. Also the former present much higher compressive and tensile strengths while keeping the same microproperties for both specimen types. That forced us to perform another calibration process of microparameters for this form of particle packing.

We obtained the first results with the coupled FLAC-PFC2D scheme using an hexagonal PFC2D inclusion 260 m length by 8.5 m height embedded in the Beaumont dolomite layer. The development of microcracking is in accordance with the predicted tensile zone obtained with a fully continuum approach. On the other hand, the sudden rise up of the microcracks number can be interpreted as the damage initiation within the competent beds. A conclusive statement about the failure mechanisms could be achieved when information about real microseismic event locations will be available.

Let us notice that this work is a preliminary one in such complex geomechanical modelling field. In addition, the presented results remain specific to several simplified assumptions taken compared to reality; in particular for the cavity geometry and its evolution, but also for the geomechanical properties of certain materials at the immediate roof and in the overburden, which influence the behaviour of the competent layers. This research opens many new perspectives in modelling the induced seismicity of the study area. The following steps should be: (a) after installation of the initial stress state in the PFC2D inclusion by assuming a more realistic particle packing (random packing), to perform a sensitive analysis on the microcrack threshold, and evaluate its effect on the mechanism of microcrack growth, while computing the associated AE; (b) to compare these numerical modelling results with those of the measurements recorded in situ after restarting of salt exploitation in a few months.


Cauvin Maxime
email: maxime.cauvin@mines.inpl-nancy.fr

Verdel Thierry
email: thierry.verdel@mines.inpl-nancy.fr

Salmon Romuald
email: romuald.salmon@ineris.fr

MODELLING UNCERTAINTIES IN PILLAR STABILITY ANALYSIS

Keywords: Uncertainties, Risk analysis, Monte Carlo simulations, Pillar stability Incertitudes, Analyse de risque, Simulations de Monte

Many countries are now facing problems related to their past mining activities. One of the greatest problems they have to deal with concerns the potential surface instability. In areas where a bord-and-pillar extraction method was used, deterministic methodologies are generally used to assess the risk of surface collapses. However those methodologies suffer from not being able to take into account all the uncertainties existing in any risk analysis. Through the practical example of the assessment of a single pillar stability in a very simple mining layout, this paper introduces a logical framework that can be used to incorporate the different kinds of uncertainties related to data, models as well as to specific expert's choices in the risk analysis process. Practical recommendations and efficient tools are also provided to help engineers and experts in their daily work.

Introduction

Many countries are facing problems related to abandoned underground mines. In France, iron ore has been industrially mined all around the Lorraine Region (North-Eastern France) from the middle of 19 th to the end of the 20 th century. The most common mining method was pillar extraction but a bord-and-pillar partial extraction method was sometimes used under sensitive zones to protect houses and surface infrastructures from subsidence during the exploitation. In those areas, some collapses already occurred, either during the mining operation or a long time after it. Even though numerous studies have already been carried out, researchers still go on developing and improving their works in order to understand the reasons of such collapses.

Today, one of the greatest issues of the work of state organisations having to manage the consequences of those past mining activities is to deal with the problem of potential surface instability. In the Lorraine region, INERIS and GEODERIS developed a multi-criteria analysis in order to realise a risk map covering all the iron ore mining areas [START_REF] Merad | Use of a multi-criteria decision-aids for risk zoning and management of large area subjected to mining-induced hazards[END_REF]. In The Netherlands, a deterministic method was applied successfully in the recent past to investigate the necessity of underground support measures to protect the surface from collapsing above shallow abandoned bord-and-pillar limestone workings [START_REF] Bekendam | Stability and subsidence assessment over shallow abandoned room and pillar limestone mines[END_REF]. Even though those methods are perfectly efficient in their field of competence, they suffer from not being able to take into account all the uncertainties inherent to geotechnical problems.

In the past decades, literature about the use of probabilistic methods in geotechnical engineering became more and more comprehensive [START_REF] Einstein | Risk and risk analysis in rock engineers[END_REF][START_REF] Baecher | Reliability and Statistics in Geotechnical Engineering[END_REF]. For the purposes of any risk analysis, one may distinguish between four classes of uncertainty ranging from a very general to a very specific degree. Uncertainty can thus be attributed to (1) the scientific, economical and politic context of the study; (2) the expertise laying on deterministic human choices; (3) the use of models; and (4) the randomness and/or the lack of knowledge on data (figure 1). It can be shown from figure 1 that these four classes of uncertainty are not totally independent and that dealing with the uncertainty associated to a specific class may have impact on uncertainties of lower classes. However, the way of treating uncertainties is different between all the categories. In the higher classes, uncertainties are very difficult to deal within a quantitative way while tools exist to help considering lower classes uncertainties. In geotechnical engineering, most of the probabilistic studies that have been undertaken concern slope stability analysis. This geotechnical field indeed offers a really efficient framework for the incorporation of uncertainty into slope design [START_REF] El-Ramly | Probabilistic stability analysis of a tailings dyke on presheared clay-shale[END_REF] because physical phenomena are relatively well known and modelled. Nevertheless probabilistic studies are relatively rare in the mining field. Efforts have started to be made for the last couple of years in New-Zealand [START_REF] Richards | Waihi Underground Workings Stage II Investigations (Risk Assessment & Mitigation)[END_REF] or in Great Britain [START_REF] Swift | Stability problems associated with an abandoned ironstone mine[END_REF] after surface collapse events in order to determine where further such collapses are likely to occur in the future. Even though those studies give really interesting results and make it possible to quantify the risk of future surface collapses, they do not take into account all the different kinds of uncertainties.

The study presented herein focus on the practical assessment of a single pillar stability in a very simple mining layout. It provides a methodology to incorporate different kind of uncertainties related to data, models or specific expert's choices. Different engineered recommendations will thus be given.

Classical deterministic approach

Bord-and-pillar method of mining is commonly used to extract ore from horizontally bedded sedimentary strata. Pillars are often left intact to provide permanent support to the undermined roof. In such a context of exploitation, economical considerations and safety assessment are the two most important issues at stake. Over the past few decades, several works about pillars have been carried out by a number of researchers and engineers either to determine optimal design pillar sizes [START_REF] Salamon | A study of the strength of coal pillars[END_REF][START_REF] Van Der Merwe | Linear Coal Pillar Strength Formula for South African Coal[END_REF] or to assess long-term pillar stability [START_REF] Bekendam | Stability and subsidence assessment over shallow abandoned room and pillar limestone mines[END_REF].

The conventional methods used to assess the stability of mine pillars are based on deterministic (empiric or analytic) approaches. They usually adopt the safety factor (SF) as an indicator of the stability of the pillar. This safety factor is defined as the ratio of the pillar strength (R) over the mean vertical stress acting on pillar (S). Theoretically, a SF value greater than 1 thus means the system is stable while a SF value lower than 1 means it is unstable. In practice, a threshold value higher than 1 is generally used at the design level to incorporate uncertainties.

Several methods can be used to determine both of S and R. S is generally calculated using the Tributary Area Theory, which considers the total overburden load directly over the pillar and the portion of the galleries at its perimeter:

lL ) L b )(L l b (l γD S + + = , (1)
where l and L are the pillar width and length, b l and b L the bord width and length, γ the average volumic weight of the overburden and D, the mining depth.

Pillar strength can be calculated by laboratory measurements on samples or by back analysis of past collapses. In the Lorraine region, several field investigations have been made after the historical surface collapses events and a back-analysis of the data lead experts to evaluate the long-term strength of the pillars to be 7.5 MPa.

In the present approach, the mining layout being studied is very simple. It consists in a grid of 3 by 3 rectangular pillars that allows to take into account the mining environment of the central pillar (figure 2). Dimensions of pillars and characteristics of the mining area have been extracted from a huge database containing information about all the Lorraine mines. The exploitation panel selected for the purpose of the study is between 130 and 150 m deep and 500 m wide. Pillars have a 12 by 8 m rectangular shape, roads are 4 m wide and the seam is 4 m high. The overburden is assumed to have a mean volumic weight of 25 kN.m -3 . Using the tributary area theory and a constant strength value of 7.5 MPa, the pillar is assessed to be stable with a safety factor equal to 1.09.

Taking into account the uncertainty on input parameters

Considering figure 1, uncertainty on input data can be linked to the spatial and temporal natural variability of mechanical and geological parameters as well as to the lack of knowledge making it difficult to assess exactly the parameter values used in models.

In our case of interest, both kinds of input data uncertainty can be found. 'Physical' parameters, as γ, are affected by natural variability. They are indeed geological averaged values obtained from sample measurements. On the other side, epistemic uncertainty is predominant for parameters as the long-term strength of the pillar (R) that has been evaluated from a statistical back-analysis of several surface collapses occurred in the iron ore basin.

The determination of the nature of the uncertainty is far more difficult for the geometrical parameters such as pillar dimensions or mining depth. These are usually determined from mining paper maps archived by mining companies or regional authorities, and which are generally at a 1:5000 scale. At such a scale, a pencil line thickness on the map is equivalent to more than 1m in the mine. For classical Lorraine pillars, such a map inaccuracy may lead to a 10 % error on the values used for pillar dimensions and obviously to significant errors in the risk analysis process. It has been considered in this study that uncertainty due to the map inaccuracy belongs to 'model uncertainties' as maps are models aiming to help engineers to represent the reality.

In geotechnical engineering, a very efficient way of taking into account uncertainty on input parameters is to define statistical distribution functions for each of the parameters [START_REF] Kim | Probabilistic Approaches to Estimating Variation in the Mechanical Properties of Rock Masses[END_REF]. Unfortunately, available information is frequently poor and test data are often inefficient to characterize the statistical moments required in the distribution functions. When test data are unavailable, those parameters may be estimated from literature or judgement and experience of experts (USACE, 1999).

Once distribution functions are defined, Monte Carlo simulations may be used to propagate uncertainties in the deterministic model. The principle of such simulations is to take one random value from the distribution of each input parameter and to calculate the output value of the function. Each calculation gives one outcome. A great number of simulations allows to obtain a frequency distribution of the output values. Mean value, standard deviation or other statistical moments of the studied parameter may thus be estimated. The greater the number of simulations, the more accurate the computed statistical moments.

For this study, a simulation method has been built allowing all the parameters to be defined by a mean value and a standard deviation. This method has been computed using Mathematica software.

For the construction of the mining layout, each position of the corners of the nine pillars is randomly distributed around a mean value. In order to respect mine design reality, the assumption was made here that the roads are cut in a straight line. Map inaccuracy is supposed to lead to a 1 m error on pillar dimensions. Figure 3 presents different layouts that can be obtained from random simulations. For illustrating the difference between them, safety factors (SF) and extraction ratios (τ) were calculated for a 140 m thick overburden. In order to describe statistically other input parameters, several assumptions had been made. All parameters had been characterized by a mean and a maximal error. Normal distributions were chosen to fit all the parameters. Such an assumption may be valid for parameters that have been determined from back-analysis of several past events (as R) or that are considered to be traditional mean values of sedimentary rock properties (as γ). It is more opened to criticism for parameters, as h or D, that are described in the available database by a minimal and maximal value. However, it will be shown later that the choice of a specific distribution function has only a weak influence on the result (see table 2). The maximal errors chosen for the parameters are considered to be the boundaries of the 95 % confidence interval of those normal distributions. When no information was available, maximal errors were roughly assumed to be 10 % of the mean value. Table 1 summarizes the statistical values chosen for each of the input parameters. Using the previously described method, the statistical analysis of a large number of computed simulations made it possible to evaluate the mean value and the standard deviation of the safety factor. Figure 4 presents the obtained distribution function. From a practical viewpoint, it would be of interest to estimate p f , the probability of 'design failure', defined as the probability that the computed safety factor is less than the target value of stability, i.e. 1:

) 1 SF ( ob Pr p f < = . ( 2
)
p f is in this case equal to 23.5 %. While using of Monte Carlo method, an important parameter to take into account is the number of computed simulations. The greater the number of simulations, the more accurate the computed statistical moments but the longer the computation. In this study, a parametrical study has been carried out in order to investigate this parameter. Figure 5 presents for different numbers of simulations the ranges of variation of computed means of safety factors, as well as times to compute them. It can be shown that, for 10000 simulations, the computation time is still reasonable and the range of variation of the computed safety factor is less than 1 %.

Taking into account the uncertainty related to the choice of a specific model

As presented earlier, uncertainty in risk analysis can also be attributed to some of the expert's choices. Assessing a risk generally requires using specific methodologies and mathematical models. For example, in our current example of pillar stability assessment, the choice has been made to compare the strength of the pillar to the stress acting on it. Two models have thus been used to estimate both of those parameters. Actually, models are no more than abstractions of the state of nature and no matter how sophisticated they are, they are unable to capture the nature entirely. The most effective approach to estimate the previously named 'model uncertainty' would be to rely on field observations and use databases to do back-analyses. However, in geotechnical studies, it is usually very difficult to have historical observations to compare with results of model predictions and most of the times, only the uncertainty related to the choice of specific models (classified as 'expert uncertainty') is taken into account.

In a study by Husein [START_REF] Malkawi | Uncertainty and reliability analysis applied to slope stability[END_REF], uncertainty associated with the use of different slope stability models was addressed by evaluating the relative performance of some simplified methods compared to the most rigorous and accurate one. In this study, it has been chosen to consider models as input parameters. Such a choice allows taking into account both the existence of several mathematical models and the inability of the analyst to identify the best one.

In our work previously described, two mathematical models have been used: the tributary area method has been chosen to calculate the stress acting on the central pillar of the mining layout and a constant value model has been used to evaluate the strength of a pillar. Nevertheless, even if both of the models are really convenient because simple, they also have disadvantages.

The two main drawbacks of the tributary area method are that (1) it makes the assumption that each pillar is carrying the load of a vertical rock column over it, this may only be valid for central pillars from large horizontally stratified mining environments, and (2) it does not take into account the nature of the overburden. Even though numerical simulations such as finite elements methods are today the most common methods used to compute stresses acting on pillars, they can only be used for specific sites where geology and mining geometry are known. Nevertheless, for the purpose of this study, whose aim is to be as general as possible, it has been chosen to use results from Coates' work about stresses acting on pillars [START_REF] Coates | Principes de la mécanique des roches[END_REF]. His approach makes it possible to take into account the extent of the mined area, the stress component parallel to the seam, the relative deformation properties of pillar, roof and floor rocks, and the positions of the pillars in the mining zone. The general solution for the mean pillar stress becomes:

( ) ( ) ( ) ( ) ( ( ) ) ⎟ ⎟ ⎟ ⎟ ⎟ ⎠ ⎞ ⎜ ⎜ ⎜ ⎜ ⎜ ⎝ ⎛ ν - ν - τ + ⎟ ⎠ ⎞ ⎜ ⎝ ⎛ + τ - + ν - ν - ν - ν - - τ + γ = w w p w p w 0 p p w w 0 2 1 2 1 A l 2 N 1 1 1 2 E E A h E E A h K 1 1 2 1 A h K 2 1 D S (3) with ) b l )( b L ( Ll 1 + + - = τ , ( 4
)
where τ is the extraction ratio, K 0 the coefficient of geostatic stress, A the extend of the mining area, N the number of pillars, E p , E w and ν p , ν w the Young modulus and Poisson ratios of pillar and wall (roof and floor) materials.

For the purpose of the study, the assumption had been made that only the parameters that are present both in equations ( 1) and ( 3) are aleatory. The impacts of those equations on the result may thus really be compared. Table 1 presents the values chosen for all the parameters. Some of these had been extracted from local databases and some are roughly estimated.

In our case of study, the use of the tributary area theory leads to a stress value of 6.87 MPa while Coates' formula gives a value of 6.66 MPa. It confirms the general belief that the tributary area theory is rather a conservative approach.

Concerning the strength of a pillar, it is now common cause that it depends on the strength of the material of which it is composed as well as its dimensions, and more specially its width-to-height ratio [START_REF] Salamon | A study of the strength of coal pillars[END_REF][START_REF] Van Der Merwe | Linear Coal Pillar Strength Formula for South African Coal[END_REF]. Using Van der Merwe's method on a very extensive database of failed and stable Lorraine iron ore pillars, it was found that the strength of those pillars can be expressed as follows [START_REF] Cauvin | Elaboration d'une méthode permettant de distinguer risque d'effondrement brutal et risque d'effondrement progressif des terrains suite à une exploitation du sous-sol[END_REF]:

h w k R 2 = (5) with l L Ll 2 w + = . ( 6
)
Here, w is the equivalent width that has to be used for rectangular pillars and k a constant value that can be related to the strength of the pillar material. The determination of k lays on the assumption that, on average, pillar failure occurs for a safety factor equal to 1. Figure 6 shows the comparison between the two different models (constant-value versus equation ( 5)). For a w/h-ratio higher than 2.9, R 2 is bigger than 7.5 MPa.

Those different physical models have been implemented in the simulation method. Basically, like for input parameters, one random model is chosen among the two that exist to evaluate the strength of the pillar and the stress acting on it. The computation is then done following the same procedure as described in part 3. A degree of confidence can be chosen for models allowing taking into account decisions of experts as well as historical considerations. For example, if in a panel of 5 experts, 3 think that it is better to use a constant-strength value to express the strength of the pillar, the strength model used in the computation have a probability of 3/5 (0.6) to be this one. In one other way, in the Lorraine region the tributary area theory was originally used by miners to design the underground workings. It thus seems logical to preferably use this formula to assess the risk of pillar failure. Figure 7 shows the results obtained by the simulation method, assuming that the input parameters are perfectly known. It has been chosen here to assign a confidence degree of 75 % for the tributary area theory while the models for the pillar strength are equally "trusted". It can then be shown in this case that one on the four combinations gives a safety factor less than 1. This Safety Factor corresponds to the ratio between R 2 (equation 5) and S (equation 1).

Figure 8 presents the distribution function obtained for the safety factor when uncertainties on input parameters and on models are simultaneously taken into account. In comparison with the previous results, the probability of 'design failure' increases from 23.5 % to 34 %.

Taking into account the method used to "propagate" uncertainties

Monte Carlo simulations have been chosen as the method used to propagate uncertainty on input parameters into deterministic models. We are fortunately here in a case within which it is possible to do so. Problems with Monte Carlo simulations are that they may require huge time of computation.

However others methods exist to "transfer" the uncertainty concerning the input parameters on the output of a model. The Point Estimate Method (PEM) is one of them. It is known for its simplicity and is commonly used in geotechnical studies. One of the greatest advantages of such a method is that it is easily computable in numerical simulation processes. In PEM, each random variable is defined by two numbers: the mean value plus or minus the standard deviations of the variable.

Using those points, it allows to obtain accurate approximations of the mean value and the variance of unknown output variables. For more details about the method, a detailed algorithm is presented in [START_REF] Isaksson | Model for Estimation of Time and Cost Based on Risk Evaluation Applied on Tunnel Projects[END_REF]. The main drawback of PEM is that it does not provide any information about the shape of the distribution of the output parameter. Moreover, using complex functions or numerous correlated variables may lead to an inaccuracy of the results.

In our case of interest, Monte Carlo simulations allowed to obtain the distribution function of the pillar safety factor (figure 4). The mean value and the standard deviation of this parameter can thus be estimated to be respectively 1.095 and 0.126 after 10000 simulations. Using PEM, they have been found to be respectively 1.104 and 0.163.

Analysis of variance

Engineers and experts clearly ask today for efficient tools that can help them to solve the tough dilemma they face. They indeed have to give quick and accurate answers to problems full of uncertainties. One of their greatest needs is thus to reduce the uncertainty on the results they give (in this study, the standard deviation of the calculated safety factor is about 0.14) in a rapid and cheap way. For time and economic reasons, it is obvious that engineers cannot study all the sources of uncertainty but can only try to focus on some of them. The study presented herein provides a tool that makes it possible to identify the parameters whose uncertainty has the greatest impact on the result.

Table 2 presents the relative influences of five different sources of uncertainty on the total variance of the computed safety factor. The map inaccuracy and the use of human-made database, expertmade hypothesis or models (as back-analysis) to estimate parameter values may have impacts on the input parameters. But, as presented in part 4, the choice of specific models to deal with a problem may also have impacts on the result. In order to obtain Table 2, each source of uncertainty has been investigated separately. Normal and uniform distribution have been used in the Monte Carlo simulations to fit the input parameters. It can be shown that the choice of a specific distribution function has only a weak influence on the result. It can also be highlighted that in our case of interest, almost half of the total variance of the computed safety factor is due to the map inaccuracy. It has also to be specified here that the differences between the relative influences of D, γ and R on the total variance can be totally explained by the differences between the maximal error on mean ratios as the rule they play in the calculation of the safety factor is exactly the same (see equation 1).

Conclusions and practical recommendations

Table 1 introduced a possible scheme to classify the different categories of uncertainty that can be encountered in a risk analysis study. This paper provides a logical framework that can be used to reduce several kinds of uncertainties. Defining statistical distribution functions for input parameters may help engineers to deal with the problem of 'data uncertainties' inherent to the study of a natural system. It can also be pointed out here that the definition of the characteristics of those distribution functions requires the expert to make some choices (evaluation of mean values, maximal errors, …), so that the 'expertise uncertainty' is also involved. Choices had also to be made regarding the nature of the distribution functions fitting the parameters. In our study, normal and uniform distributions have both been used and it can be shown that such a choice has only a weak influence on the results. Table 2 illustrates that the variance of the computed safety factor is only slightly bigger when uniform distributions have been chosen for all the input parameters (in this case, p f equals 24.6 %) than when normal distributions have only been used (p f = 21.4 %).

More generally about 'expertise uncertainty', this work may guide engineers through their analysis process in giving them a methodology to follow. It thus allows to break free from the individual ability or from the deterministic choices of the expert in charge of the analysis.

It has been said in introduction that the methods currently used to assess the risk of surface collapse do not generally integrate the taking into account of all the uncertainties inherent to the study of a natural system. In fact, it is not exactly true. Most of the times, uncertainties are introduced in the evaluation process in an indirect, if not unconscious, way. This integration is basically part of the expert's analysis process and is affected by 'expertise uncertainty'. For instance, in the example we used, finding a safety factor greater than 1 means the system has to be stable. However the value of 1.092 is very close to 1, meaning that we are thus very close to the stability. Experts usually use deterministic threshold values for safety factor (1.2/1.3/1.5/3) in order to integrate uncertainties inherent to the evaluation and to be side of the safety. Generally, the more damaging the consequences of an event, the greater the chosen threshold value.

Nevertheless, using such threshold values does not allow the expert to assess 'how stable the situation is'. Taking explicitly into account uncertainties on input parameters by the use of distribution functions makes it possible for the expert to represent himself an order of magnitude of what he is talking about, e.g. in our case of interest 23.5 % of the computed safety factors are less than 1. Expressing such a quantified result now sets the problem of its interpretation. In fact, we are now talking in terms of risk and acceptance. It is now up to the stakeholders, not either to the engineers, to decide whether the situation is acceptable or not, and whether uncertainty on the result has to be reduced or not.

Figure 1 .

 1 Figure 1. Possible scheme distinguishing between the different categories of uncertainty in engineering risk analysis.

Figure 2 .

 2 Figure 2. Extract of a mining map and selection of the nine pillar mining layout.

Figure 3 .

 3 Figure 3. Three mining layouts computed from random simulations.

Figure 4 .

 4 Figure 4. Distribution of the safety factor when uncertainties exist on input parameters.

Figure 5 .

 5 Figure 5. Range of variation of results and computation times for different numbers of simulation.

Figure 6 .

 6 Figure 6. Comparison between the two different models estimating pillar strength.

Figure 7 .

 7 Figure 7. Computed safety factors, taking into account the existence of different models.

Figure 8 .

 8 Figure 8. Distribution of the safety factor when uncertainties exist both on input parameters and models.

Table 1 .

 1 Statistical values used for input parameters.

	Parameter		Symbol Mean value Maximal error Distribution law
	Pillar length	(m)	L	12	1	Uniform
	Pillar width	(m)	l	8	1	Uniform
	Bord length = bord width	(m)	b	4	1	Uniform
	Pillar height	(m)	h	4	0.5	Normal
	Mining depth Average volumic weight	(m) (MN.m -3 )	D γ	140 0.025	10 2.5e -3	Normal Normal
	Long-term strength	(MPa)	R	7.5	0.5	Normal
	k-Strength constant	(MPa)	k	4.4	(-)	Constant
	Mining span	(m)	A	500	(-)	Constant
	Number of pillars	(-)	N	31	(-)	Constant
	Coefficient of geostatic stress	(-)	K 0	0.5	(-)	Constant
	Young modulus of pillar material (GPa)	E p	13	(-)	Constant
	Young modulus of wall material (GPa) Poisson ratio of pillar material (-) Poisson ratio of wall material (-)	E w ν p ν w	5 0.33 0.33	(-) (-) (-)	Constant Constant Constant

Table 2 .

 2 Relative influence of the different sources of uncertainty on the computed safety factor.

	Variance of the
	computed safety factor


Cojean Roger
email: roger.cojean@ensmp.fr

Nicolas Franco
email: 1.nicolas.franco@charbonnagesdefrance.fr

Lazarewicz Jean-Claude
email: jean-claude.lazarewicz@charbonnagesdefrance.fr

Blachere Agnès

Lefort Dorothée

Claudia Sorgi
email: claudia.sorgi@ineris.fr

): ANALYSIS OF IMPACTS RESULTING FROM THE CESSATION OF PUMPING AT THE CENTRAL SHAFT. SURVEY OF VARIOUS SCENARIOS RELATED TO THE WATER LEVEL OF THE PIT LAKE IN THE GRANDE DECOUVERTE

d'enseignement et de recherche français ou étrangers, des laboratoires publics ou privés.

ABSTRACT: Within the frame of the Survey related to the cessation of mine workings in Decazeville-Firmi concession, various impacts resulting from the cessation of pumping at the Central Shaft were assessed. Mainly these impacts are related to groundwater behaviour in the abandoned underground coalmines, hydrochemistry of waters discharged to the environment, ground stability concerns and coal gas emanations. This analysis allowed the choice of the most appropriate elevation of the pit lake level in the Grande Découverte, with the necessity to continue the pumping. Two main objectives were reached. The elected elevation is high enough to result in a permanent chemical stratification, which allows the pumping of the superficial waters and its discharge to natural watercourses without any treatment. The elected elevation is low enough to avoid any problem of inflow of water or ground stability at some particular places which might have been threatened by the rising of the piezometric level of the mining aquifer. Lastly, the elected elevation of the pit lake allows a quality scenery design around the pit lake.

Introduction

The paper is devoted to the analysis of the post-mining context at Decazeville-Firmi concession (Aveyron, France). Complex technical conditions resulting from the past coalmine workings, discussions related to the assessment of various potential impacts and the choice of the most appropriate level of the pit lake in the Grande Découverte are at the origin of a long procedure which took place during four years, from 2001 to 2004. During this time, dialogues occurred between Charbonnages de France, its engineering and technical consultants (mainly Ecole des Mines de Paris -Centre de Géologie de l'Ingénieur, CESAME, INERIS), the Midi-Pyrénées DRIRE (Regional Direction for Industry and Research), the local and regional authorities (Villefranche-de-Rouergue Sous-Préfecture and Rodez Préfecture) and the elected representatives of the concerned communes. After several investigations in the fields of hydrogeology, engineering geology and geomechanics, and presentation of the obtained results to the various concerned authorities, a general technical and consensual agreement was obtained. In a brief way, the paper describes the main steps of the corresponding scientific or technical investigations which allowed to address and solve the impact problems (Cesame, 2000(Cesame, , 2004;;Charbonnages de France, 2001;[START_REF] Cojean | Compte-rendu d'expertise[END_REF][START_REF] Cojean | Note complémentaire au Compte-rendu d'expertise relatif au dossier d'arrêt des travaux miniers dans la concession de Decazeville-Firmi (Aveyron)[END_REF][START_REF] Cojean | Note complémentaire au Compte-rendu d'expertise relatif au dossier d'arrêt des travaux miniers dans la concession de Decazeville-Firmi (Aveyron)[END_REF][START_REF] Ineris | Bassin houiller d'Aubin-Decazeville[END_REF][START_REF] Ineris | Bassin houiller d'Aubin-Decazeville[END_REF][START_REF] Lefort | Etude hydrogéologique du bassin houiller de Decazeville (Aveyron). Thèse professionnelle[END_REF].

Regional geological setting

According to geological investigations mainly due to Vetter P. ( 1968) and later research works resulting from the coal mining works, the litho-stratigraphic sequence of the coal basin, as well as the tectonic deformation phases are well known (Figures 1 and2). The sedimentary sequence consists of several cycles or "assises" dated from Middle Stephanian series (Carboniferous period). In a stratigraphical order, the cycles are arranged as follows: Brayes cycle, Volcanic complex, Auzits cycle, Banel cycle, Campagnac cycle and Bourran cycle. Each cycle is made up of conglomerates, sandstones, shales and coal layers. The major tectonic phase occurred at the end of the Paleozoic era and resulted in a general compression of the coal basin, with the folding and faulting of geological strata (Figures 2 and3).

Past coalmine workings

Decazeville coal mining basin experimented three major mining periods. From 1360 to 1744, numerous private exploitations developed. Then, from 1744 to 1946, the period of large mining companies took place. In 1946, nationalization of mining companies and creation of Charbonnages de France occurred. It was a period of intense mining activities, with peak productions, then decline due to the international economic coal context and progressive lowering of the regional coal resources (Charbonnages de France, 2001). Underground mining works were performed according to a large variety of mining methods (Figure 4). Some of them were partial exploitations such as the method of chambers and pillars that left important residual voids at the end of the exploitations in case of abandonment. Other ones were total exploitations with an integral void treatment such as the stoping methods by successive horizontal slices: the cut and fill method that is working upwards or the undercut and fill method that is working downwards.

During the recent period of large opencast exploitations and especially the progressive deepening of the Grande Découverte, as well as other open pits in the surroundings, old mining works were partly intersected. Field analyses allowed to investigate the old backfilling operations as well as the mechanical behaviour of grounds after mining subsidence.

Post-mining conditions

The post-mining context was analysed in a detailed way, considering the following causes for concern and corresponding impacts on the environment and constructions [START_REF] Feuga | Quelques aspects de l'après-mine[END_REF][START_REF] Fabriol | La gestion de l'eau des mines en phase post-extractive[END_REF]:

̵ Type and behaviour of the mining aquifer;

̵ Ground stability concerns;

̵ Impacts of various scenarios related to the cessation of pumping and the water level in the pit lake; ̵ Risks resulting from mining gases, surface mining facilities, waste dumps and tips, specific equipments.

Only the first three items are presented in this paper.

The mining aquifer

The question of underground waters was investigated considering the different types of aquifers as well as their properties by reference to the problem of discharges and the question of water pressures and ground stability (Cesame, 2000;[START_REF] Cojean | Compte-rendu d'expertise[END_REF][START_REF] Younger | Mine water[END_REF][START_REF] Schmitt | Qualité des eaux après fermeture des mines : remplissage initial, évolution transitoire, stabilisation à long terme et gestion environnementale[END_REF]. The Carboniferous rock masses can be considered as made up of three different hydrogeological units that were more or less intensively modified by the mining works and associated mining subsidence. The Table 1 summarizes the characteristics of these units.

As a matter of fact, the water-bearing formation can be considered as mainly determined by the old coalmine workings, generally backfilled by permeable materials, while discharge questions and water inflows are concerned. However, for underground water pressures and stability problems, the hydraulic continuity has to be underlined as the major concern and may result from continuous water paths through low or very low permeability formations. So, in order to completely investigate these questions, very detailed consultations of old mine working drawings were performed in order to understand the connections between different mining works or to identify the presence of dams in particular cross-cut drifts.

In addition, a piezometric network was installed to get a good understanding of the underground water behaviour and improve the knowledge of water pressure distributions in some specific places.

Table 1. Hydrogeological behaviour of the mining aquifer. Storativity and transmissivity are highlighted. In addition, the exchange function related to interactions between underground water and storage rocks has to be considered for physico-chemical, hydro-biological and thermal processes. It controls water quality.

Mining aquifer Storativity function Transmissivity function

Carboniferous rock masses directly related to the old coalmine workings. Very high permeability. Very high porosity.

Very significant role Very significant role

Carboniferous rock masses altered by the old coalmine workings and subsidence phenomena.

Medium to low fracture permeability. Medium to low porosity.

Moderately significant role Poor role

Carboniferous rock masses unaltered by the old coalmine workings.

Low (subsurface zones) or very low permeability. Low (subsurface zones) or very low porosity.

Very poor role Very poor role

As an illustration, conclusions resulting from the analysis of data related to two piezometers are presented in Figure 6. In the two cases mining connection was existing, but it can be seen that the hydraulic connection does not have the same reaction while the hydraulic gradient is lowering due to the rising of the water level in the Central mineshaft (Puits Central) which is used as a pumping shaft. Taking advantage of these different investigations a synthetic map was provided (Figure 7) which illustrates the different hydrogeological units and connections with the Central mineshaft.

The ground stability concerns

Ground stability questions were investigated by reference to potential mining subsidence, collapse of abandoned rooms and pillars exploitations and sinkhole processes [START_REF] Ineris | Bassin houiller d'Aubin-Decazeville[END_REF][START_REF] Piguet | Affaissements et déformations au-dessus des exploitations minières : mécanismes et évolution dans le temps[END_REF][START_REF] Cojean | Compte-rendu d'expertise[END_REF][START_REF] Wojtkowiak | Maîtrise de l'impact hydraulique, sur le bâti de surface, de l'ennoyage de mines souterraines arrêtées du bassin ferrifère lorrain[END_REF][START_REF] Heib | Paramètres d'affaissement pour la hiérarchisation des zones à risque dans le bassin ferrifère lorrain[END_REF].

The areas which formerly were affected by mining subsidence were delimited. Presently, no residual subsidence can be expected in these places. However, reactivation of mining subsidence as a response to the rising of the general water level in the mining aquifer had to be considered. Different observations in the world show that some low uplift process can be expected. It has to be noted down that these processes occur at a regional scale. On the contrary, at a local scale, the rising of underground water level can result in local reactivation of settlement due to the washing out of infilling materials in sloping galleries for shallow exploitations.

The collapse of partial exploitations such as rooms and pillars is a concern for exploitations with a depth lower than 50 meters. The present stability conditions were investigated using the data provided by colliery archives. Long term stability conditions were analysed, taking into account the lowering of strength of pillars due to the rising of water levels in the exploitations. It was necessary to carry out additional boreholes in some places in order to clarify particular situations where safety problems were at stake. The sinkhole hazard was also investigated considering possible scenarios of cavities at the crown of an underground mine working turning into successive bell holes reaching to the surface. For this purpose, additional boreholes have allowed to precise the risks associated to these processes and gave some essential data for stability calculations. Lastly, ground stability questions were also analysed for underground iron mines, connected to the coalmine workings, which would be concerned by the rising of the regional piezometric level as a response to the cessation of pumping at the Central mineshaft.

As it can be noted down, any ground stability problem is closely connected to underground waters and corresponding water pressures. That is the reason why it is necessary to analyse the roles of mining aquifers in post-mining conditions by respect to discharge and water inflow questions, as well as by respect to piezometric levels and water pressures.

Election of the most appropriate water level in the pit lake

At the beginning of the Survey related to the cessation of mine workings in Decazeville-Firmi concession, comparative scenarios were analysed. At the beginning of the investigations, the scenario of a lake level at the 211 m elevation was considered. It was corresponding to a total cessation of pumping and made it possible to have a gravitational flow of water towards the surface watercourses, but consequences concerning ground stability concerns and water inflow in some places were not completely solved. Then a scenario of a lake level at the 170 m elevation was analysed. It allowed to cancel any ground stability problem, but long term pumping had to be taken into account as well as water treatment.

Finally an intermediate solution appeared to be appropriate. It was established according to the following thinking. The option had to make it possible to discharge water towards the surface watercourses without any treatment. But potential ground stability problems had to be limited and residual cases to be solved, especially for two cases related to Vialarels area in the town of Decazeville and Combes area in a nearby village. A consequence was that pumping had to continue. Some world expertise concerning the behaviour of pit lakes is available as well as French expertise concerning coalmining works, especially at Blanzy (Saône et Loire). Pit lake evolution is complex and depends on geometrical parameters related to the watershed, hydro-chemical parameters concerning the lake supply from surface and underground flows, characteristics of the connection between the lake and the mining aquifer and climatic parameters. The seasonal thermal stratification is well known for lakes in general. But for mine lakes, a chemical stratification can also take place.

A permanent confined and mineralized body of water appears at the bottom of the lake and does not participate to the seasonal thermal convection processes. This zone appears below a narrow transitional zone called the chemocline, the corresponding lakes being described as meromictic lakes. According to the expertise coming from Blanzy [START_REF] Paquette | Les lacs miniers des exploitations à ciel ouvert de charbon du massif central : qualité et évolution[END_REF], three main factors seem to control this chemical stratification with threshold values: ̵ The maximum depth of the lake has to be larger than 25 m;

̵ The mean depth (ratio volume/surface: V/S) has to be larger than 15 m;

̵ The hollow index (1000×V/S 3/2 ) has to be larger than 50.

In addition, two additional factors may play a significant role: the ratio watershed surface/lake surface (BV/S) and the age of the lake.

As a consequence, in a first step, various scenarios concerning the water elevation in the pit lake were studied and led to the results presented in Table 3.

According to the accepted threshold values, it appears that the most appropriate water elevation is 195 m, taking into account a small safety margin (Figure 10). Then it will be possible to pump a superficial water slice within a specific structure designed to pump water with no disturbance of the chemical stratification of the lake. In a second step, the consequences related to ground stability problems resulting from the 195 m elevation were considered.

The Vialarels area in Decazeville was closely analysed. This area consists in three different levels of galleries (Figure 8-a), the highest voids being locally at a 6 meters depth. These mining works were completed in 1916, and different small mining damage occurred during and after mining: small settlements, gas emanations, minefires. Infilling and reinforcement works were carried out. The hydraulic connection between Vialarels area and the mining aquifer connected to the Central mineshaft (where pumping was taking place) was specifically analysed using maps and additional piezometers which were progressively set up. The hydraulic connection appeared to be very poor (Figure 6-b), especially due to the presence of a permanent dam erected across the single cross-cut gallery connecting Vialarels area and the rest of the mine. In fact the water table in this area is mainly dependant on the recharge corresponding to direct meteoritic waters, but there is also a small flow which leaks (probably through the dam) to the mining aquifer connected to the Central mineshaft. This leak can probably be influenced by the level in the pit lake if it is upper than 195 m. That is the reason why a pumping station will be designed to pump water out in case of heavy rainstorm on this area. Finally additional analyses were performed about some other places, limited in extension, were similar impacts were identified. As a consequence of these investigations in the field of mining groundwater and ground stability, the 195 m elevation was accepted as the best solution for the water level of the pit lake (Figure 10). In addition this level allows a quality scenery design around the pit lake.

Conclusion

As a conclusion, different major subjects that were investigated in the case of Decazeville-Firmi concession can be underlined.

The post-mining contexts are always complex and need integrated investigations mainly in the field of mining groundwaters, ground stability and gas emanations. The questions of underground waters have to be studied considering the discharge as well as the water pressure problems.

The monitoring of the mining aquifer by piezometric networks is necessary in order to improve the understanding of these complex groundwater systems.

The reference to colliery archives is essential to obtain a good knowledge of past mining works.

Unfortunately the corresponding data are sometimes partial or difficult to be interpreted. So it is necessary to consider that boreholes investigations are sometimes necessary in order to accurately describe the subsoil conditions.

In the field of pit lakes, expertise begins to be accumulated but it is necessary to be careful and to implement monitoring devices in order to follow the chemical stratification process, as it is done in the case of Decazeville.

  The post-minig context at Decazeville-Firmi concession (Aveyron, France) : analysis of impacts resulting from the cessation of pumping at the central shaft. Survey of various scenarios related to the water level of the Pit lake in the Grande Decouverte Roger Cojean, Nicolas Franco, Jean-Claude Lazarewicz, Agnès Blachere, Dorothée Lefort, Claudia Sorgi To cite this version: Roger Cojean, Nicolas Franco, Jean-Claude Lazarewicz, Agnès Blachere, Dorothée Lefort, et al.. The post-minig context at Decazeville-Firmi concession (Aveyron, France) : analysis of impacts resulting from the cessation of pumping at the central shaft. Survey of various scenarios related to the water level of the Pit lake in the Grande Decouverte. Symposium Post mining 2005, Nov 2005, Nancy, France. pp.NC. ineris-00972519

Figure 1 .

 1 Figure 1. Geological cross section of Decazeville coal basin. The Stephanian series is represented in the central part of the cross section. La Découverte is concerned by the Bourran (h 5VI ) cycle. Old underground coal works are mainly related to the Bourran (h 5VI ) and Campagnac (h 5V ) cycles. (after BRGM. 2001. Decazeville geological map, 1/50000).

Figure 2 .

 2 Figure 2. Elisabeth open pit, Firmi. Photo dated 1896. Saint-Eugène anticline, Campagnac cycle. (in Vetter P., 1968).

Figure 3 .

 3 Figure 3. NE-SW geological cross section. Lassale (Grande Découverte) and Sérons open pits. Bourran cycle. (after Vetter P., 1968).

Figure 4 .

 4 Figure 4. Detailed old mining maps. a) The Bourran map is related to mine workings dated before 1893, with detailed drawings concerning backfilled horizontal slices during the period 1893-1897. b) The Paleyrets map is related to mine workings dated from 1830 to 1947, with detailed drawings concerning exploration galleries, abandoned coal chambers and pillars. In addition, the old Tramont abandoned iron chambers and pillars can be noted down. (after Charbonnages de France archives).

  a) Old gallery backfilled by slime materials b) Old timbered and backfilled mining works

Figure 5 .

 5 Figure 5. Intersections of old mine workings during the deepening of the Grande Découverte.

Figure 6 .

 6 Figure 6. Hydrogeological connections between several coalmine workings and the Central mineshaft (Puits central).

Figure 7 .

 7 Figure 7. Synthetic map showing different mining aquifers connected to the Central mineshaft (Bourran-Lassalle, Paleyrets and Buscalie sectors) and other mining aquifers with no connection to the Central mineshaft (Banel, Placards, Martinie-Moulin and Firmi sectors). (after Cesame, 2000).

  a) Mining works at Vialarels area. Three different levels of galleries can be observed. b) Rooms and pillars iron mines at Combes area. These exploitations are connected to the Central mineshaft thanks to coal galleries.

Figure 8 .

 8 Figure 8. Mining works at Vialarels area and Combes area.Also the Combes area was analysed in a detailed way. In this place, old underground iron mines are connected to the coal mines and the Central mineshaft. The three dimensional characteristics and location of the chambers and pillars iron mines were studied thanks to the colliery archive documents (Figure8-b) and knowledge related to the mining aquifer. It resulted that only a few houses were potentially concerned by modifications of stability conditions due to a variation of water level up to 195 m in the rooms and pillars of iron mines. Appropriate decisions will be taken for this particular case.

Figure 10 :

 10 Figure 10: Aerial view of the pit lake in the Grande Découverte. The Central mineshaft is on the right of the lake. In the top left corner, Decazeville town and Vialarels area. In the bottom right corner, Combes area.

Comparaison des niveaux d'eau Crucifix/Buscalie et Puits Central /Lac

		210
	Cote NGF du niveau d'eau	110 120 130 140 150 160 170 180 190 200 02/01/1999	02/01/2000	01/01/2001	01/01/2002	01/01/2003	01/01/2004 Puits Central LAC Crucifix Buscalie/Fontverhne Moy. mobile sur 2 pér. (LAC)	31/12/2004	UNITE BOURRAN-LASSALLE, secteur Crucifix-Buscalie

Table 3 .

 3 Analysis of factors controlling the chemical stratification in the pit lake.

	Water elevation	170 m	185 m	190 m	195 m	204 m	211 m
	in the pit lake
	Maximum	16	31	36	41	50	57
	depth (m)
	Volume V (m 3 )	220 000	920 000	1 330 000	1 800 000	2 800 000	3 750 000
	Surface S (m 2 )	33 000	77 000	87 000	97 000	124 000	155 000
	Mean depth V/S	6,7	11,9	15,3	18,6	22,6	24,2
	(m)
	Hollow index	36,9	42,9	51,9	59,7	64,2	61,5
	1000×V/S 3/2
	Ratio BV/S	47	20	18	16	12	10