T Geldmacher

A Studeny

Dmt Gmbh

Co

Ger M Kg

INERIS, FRA Al Heib

NUMERICAL MODELLING AND EMPIRICAL TOOL FOR SUBSIDENCE PREDICTION BY THE USE OF BACKFILL

Geomechanical-numerical methods are an established tool for planning of underground excavations worldwide. By describing the interaction between rock mass and backfill it is possible to determine the suitability of backfill methods regarding subsidence and a potential increase of the extraction ratio.

DMT and INERIS, both have longterm experience in investigating and applying geotechnical parameters for designing and refining numerical models for planning and evaluating almost any underground excavation-project by geomechanical-numerical tools combined with empirical planning methods.

The presentation describes the planning approach for numerical models in mining and illustrates the advantages of this planning tool by means of several planning cases from the mining industry worldwide. The paper also describes the development of an empirical approach visualizing the backfilling effect on the prediction of subsidence of deep mines. The estimation of the effect of backfilling may be obtained by in-situ measurements and by the numerical modeling. A basic 2D numerical model was used to estimate the factor of subsidence reduction due to the backfilling exemplary for a room and pillar mine.

INTRODUCTION Movement types relating to underground mines

The type of surface movement, the amplitude, and the extension of the movements depend on several factors: geometry and localization of the underground mine panels, depth, the nature of the overburden, the presence of water and, the mining methods: total and partial extraction methods. The prediction of the subsidence of room-and-pillar mining (partial extraction) is more complex than for e.g. longwall mining (full extraction). During the extraction the subsidence is very limited and afterwards it is difficult to determine when the mined cavities will collapse due to the degradation of pillars when mining activities have ceased.

One main function of backfilling is the generation of a counter bearing for the side and/or the roof pillars. A common method is the design of stable pillars between the rooms to guarantee the safety and to avoid rock mass fracturing and displacements (e.g. subsidence).

While implementing backfill this support effect could even be achieved with smaller pillars. The stabilization of the pillars is a result of the triaxial restraint of the pillar due to backfilling. The requirements on the material will most be little, even a loose rock in a filled case could be able to work as a counter bearing. For the stabilization of the side walls a partial backfill might be sufficient (Figure 1, middle). For a stabilization of a roof pillar a complete backfilling of the room is required (Figure 1, right). s max s max Figure 1: Room-and-Pillar with backfill A modification of the classical Room/Stope-and-Pillar method is the Room/Stope-and-Pillar method with pillar recovery. Figure 2 shows a possible layout for mining an ore body. The task of the backfill in this case is the refilling of the mined stopes to enable the pillar recovery. The requirements of the backfill depend on the dimensions of stopes and the mining sequence. Generally a cementitious backfill must be used. A low strength could be sufficient because the backfill ribs have to remain stable mainly for the period of pillar recovery.. After that a backfilling of the mined pillar area could be done with a low strength backfill.

For mining in slices (Figure 2) a stable floor or roof is required as well. In case of mining from the bottom to the top a stable floor is necessary for bearing the load of the used equipment (e.g. loader). For mining from top to bottom a stable roof is required to guarantee the safety of the mining process and the mining staff.

General backfill requirements

The backfill requirements are strongly related to the purpose of backfilling and the mining method and operation. The following parameters are essential for the characteristics of backfill: For the investigation of the influence of backfill to the mining process different planning tool are available. Basically different approaches e.g. empirical or numerical can be used. The highest accuracy can be achieved by combination of both methods.

EMPIRICAL SUBSIDENCE PREDICTION

The basic aim of all subsidence prediction methods is to produce an acceptable accurate assessment of mining subsidence (Kratzsch, 1983, Whittaker and[START_REF] Whittaker | Subsidence, Occurrence, Prediction and Control[END_REF]. Three different approaches can be applied: the empirical methods (NCB, 1975), semi-empirical (Bahuguna et al. 1993) and numerical methods (Al [START_REF] Heib | Les nouvelles méthodes de modélisation numériques et le volume d'influence des exploitations minières en conditions complexes[END_REF]. Empirical methods require feedback-data of case studies in order to predict the future subsidence for new extraction areas. This is the method of extrapolating available data to derive the future trend. Deck [START_REF] Deck | Propositions pour une méthodologie d'évaluation de la vulnérabilité du bâti[END_REF]) listed more than 30 empirical relations to predict the subsidence and the different associated parameters (strain, tilt, etc.). Empirical methods can be used to predict the maximum subsidence for different mining methods and configurations thanks to improvement, and feedback and numerical modeling, the following relationship is one of those empirical relations [START_REF] Proust | Etude sur les affaissements miniers dans le bassin du Nord et du Pas-de Calais[END_REF]:

Sm = O * f1* f2 * G
With O: the thickness of the cavity or the worked seam. For room and pillar mining the thickness of the extraction layer is replaced by an equivalent thickness equal to O'=*O. is the extraction ratio, the increasing of the extraction ratio leads to an increase of subsidence.

f1: factor which depends on the control mode of the mining area and voids, f1 varies between 0.25 (pneumatic filling) to 0.9 (total extraction and goaf). The backfilling can reduce the magnitude of subsidence.

f2: factor of depth varies from one basin to another. The value decreases with an increase of depth due to the consolidation of the caving grounds under the weight, growing with the depth of mine. Beyond the threshold, this phenomenon is counterbalanced by vault phenomena and by the fact that the conditions of appearance of sub-critical situations become increasingly frequent with great depths;

G: factor of dimensions of the exploitation, G varies according to the width and the length of the cavity compared to the depth.

As stated before in room and pillar deep underground mines, backfilling can have different functions, first to fill the void, to reduce the convergence of the strata and to increase the strength of the pillars and reduce the deformations and the subsidence. Concerning the first point, f1 factor allows the consideration of this effect. Whatever, the confining effect on pillars must be considered. The empirical equation can be modified considering the confining effect and the reducing of the convergence.

Sm = O * f1* f2 * f3 * G
Parameter f3 depends on the quality of the backfilling, the depth and the behavior of the roof. A numerical modeling was carried out to quantify the role of the confining pressure and to estimate the value of the factor f3.

Figure 3 shows a basic regular mining area with square pillars with high extraction ratio equal to 88 %. The 2D section corresponds to a deep mine at 500 m depth; the pillar section is (5 * 5 m). Instant backfilling after the extraction was supposed in the calculation.

Considering elasto-plastic behavior with Mohr-Coulomb criterion of the backfilling material, a parametric study was carried out concerning the Young modulus (50 to 500 MPa), the cohesion (0 to 1.5 MPa) and the friction angle as constant and equal to 30°. The Young modulus of the extraction layer is equal to 5000 MPa. The parameters of the backfilling can be increasing due to the increasing of the cement quantity and to the consolidation effect.

The results analyzed herein concern stress and displacement. Figure 4-a, presents the increasing of the confining pressure in the pillar and the decreasing of the convergence (on the surface) function of the increasing of Young modulus. The evaluation becomes negligible when the value of Young modulus is equal to 300 MPa: the horizontal confining stress equal to 0.6 MPa, the variation of vertical displacement (Dv) follows same tendency. Figure 4-b presents the increasing of the confining pressure and the decreasing of the convergence due to the effect of the increasing of the cohesion; the maximum effect can be obtained when the cohesion is less than 0.75 MPa. The maximum confining pressure can be reached at 1 MPa. Remember that the initial horizontal stress before the excavation is equal to 3 MPa, at the depth of the mine (500 m). The backfilling effect can be considered very important. The prediction of subsidence can be done by using empirical formula. The effect of backfilling on the stability of the pillar and the increasing confinement pressure and decreasing of surface displacement were studied using elementary numerical models. The parametric study shown with backfilling parameters of (E= 300 MPa, cohesion= 0.75 MPa, friction angle= 30°) can increase the confining pressure and decrease the convergence.

NUMERICAL MODELLING FOR DETAILED GEOTECHNICAL ANALYSIS OF BACKFILL BEHAVIOUR Numerical basics

Numerical methods are widely used for solving a wide variety of technical problems. They are mainly employed in applications where analytical and empirical processes cannot provide the required degree of accuracy because of the existence of complex boundary conditions.

As a result of this, the mining industry employs geomechanical-numerical models extensively as operational planning tools. Analytical and empirical methods have so far proved to lack the required degree of precision for describing the interaction that takes place between heterogeneous strata and various types of support elements (rockbolts, standing supports, yielding arches and backfilling) in highly deformed systems. Compared with conventional calculation methods, numerical modeling can, for example, provide a detailed analysis of the fracture and deformation status of the rock strata surrounding a mining excavation. This produces a picture of the extent of the strata fracture zone and helps dimensioning the support measures required. In order to be able to understand the complex mechanical phenomena within the rock mass, with all its inhomogeneity and separation planes, it is necessary to develop and implement various constructive models in the way of mechanical equations that can reproduce the stress-deformation behavior of the rock mass with sufficient accuracy.

DMT has used several numerical programs for the solution of different geomechanical questions for several decades. For large scaled stress redistributions the program GEDRU, developed by DMT, is used. For the most planning cases the program FLAC -2 and 3 dimensionalis applied, and in special cases other programs like PFC, UDEC or Ansys are used. The main focal point in the following is the application of the program FLAC for the solution of several geotechnical questions.

The German coal industry uses numerical models for a wide range of planning assignments, including:

• Backfilling of the goaf during longwall mining.

• Investigation of mining impact (surface subsidence, damage).

• Warranty of long term stability for underground openings like shafts, pit bottoms or main roadways.

• Investigation of the deformation behavior of temporary used underground openings like longwalls or gate roads.

• Support dimensioning.

The numerical modeling process essentially operates by reproducing all fundamental support elements (backfill, rockbolts, standing supports, roadside packs, injection material etc.) along with a wide range of boundary conditions that apply to the working of coal seams in geological deposits (geology, depth, single seam and multiple seam mining and so on).

Preparation of numerical models

The preparation of geomechanical-numerical models for control of underground openings requires a multitude of information as initial parameters. This includes information such as: To create models that are as realistic as possible, initial information with a high accuracy as well as extensive model calibration are required. To receive this information a lot of investigations have to be carried out (Figure 5): The challenge is the modification of the lab parameters for the usage in the continuummechanical program FLAC. Due to a limitation of implementing each geotechnical element (e.g. joint) in a model a limitation of these test results (e.g. rock strength) is required. Therefore different empirical approaches as well as experiences are used. An additional point is a mostly limited amount of initial parameters or their poor quality. Due to these reasons a calibration of the model is required. DMT uses 2 different ways for model calibration:

 Strata  Geology  Geotechnics  Type

 Physical modeling of underground roadways and backfill in a scale of 1:15  Underground measurements and observation

The benefit of this extensive calibration process is the output of a comprehensive picture of the complex deformation and fracture structure around underground openings. The physical modeling technique especially permits a detailed view into the rock surrounding an opening. Thus it is possible to detect fractures in the rock and the backfilled area as well as support failure (e.g. bolt failure).

Based on these considerable initial parameters and calibration capabilities DMT has successfully developed numerical modeling as a reliable planning tool which can be used worldwide. In the following examples the application of numerical models four different planning cases will be described.

Application Case

This analysis was done for a cut-and-fill operation for a coal seam deposit. Figure 6 shows the model with the boundary conditions. The rooms were headed successively. In the base model a calculation without backfill was carried out. After development of several rooms the pillars collapse completely (Figure 7). For this model several variations with different backfill rate and backfill strength were carried out. In the first variations of this basic model every room was backfilled 33% partial and 66% partial after heading parallel to heading of the next room (Figure 8 and9). Figure 10 shows the same situation as in figures 7 to 9 but with completely filled rooms. The chosen backfill compressive strength was 5 MPa for these models. The same calculations were carried out with a backfill compressive strength of 1 MPa.

The vertical displacements are an indicator for the success of the backfill regarding the avoiding the subsidence. Figure 11 shows the results of the different models. The dark blue line is the model without backfill. End of this line is the extraction of the 6th room. The green lines are the calculations with backfill strength of 1 MPa, the red/yellow lines with backfill strength of 5 MPa. Up to the excavation of the fourth room no significant vertical displacements occur due to low span width of the excavated area. No pillar collapse is indicated. In the model without backfill the pillar collapse begins after heading the next room and results in high vertical displacement. The effect of backfilling can be assessed considering the other lines. A high degree of stabilization can be achieved with full backfilled rooms as well as partial backfill with high strength. But also with a partial backfill of only 1/3 of the room and a low compressive strength of 1 MPa the vertical displacement could be reduced by 20 % compared to the unfilled room. These were the first investigations within the research project and will be amended by a multitude of variation calculations with different backfill behavior and rock mass conditions.

CONCLUSION

Based on extensive experiences over long years in underground mining and physical modeling, DMT and INERIS has accumulated a comprehensive understanding of geomechanical processes, especially regarding the interaction of rock stress, rock mass structure and support technology like backfill. This knowledge was implemented in empirical investigations and the geomechanical-numerical models. By dint of continuous calibration and optimization, a high standard of modeling is assured.

As proven by various projects this advanced modeling technique is suitable for the solution for predicting mining induced impacts such as subsidence.This confirms that numerical modeling is a versatile tool in current mine planning with increasing importance for a safe, economic and trouble-free mining operation in the future.

Figure 2 :

 2 Figure 2: Room-and-Pillar with pillar recovery and backfill

  the stowing process For the stability investigation the following mechanical backfill parameters are required: -Young Modulus, Poisson Ratio (cement-based backfill) -Bulk Modulus, Poisson Ratio (uncemented backfill) -Consolidation behavior (load-yield characteristic) -(non-cement backfill)

Figure 3 .Figure 4 .

 34 Figure 3. Numerical modeling to study the effect of the backfilling on the underground and surface behavior of room and pillar underground mine

Figure 5 :

 5 Figure 5: Required information for calibration of a numerical model

Figure 6 :

 6 Figure 6: Numerical base model with boundary conditions

Figure 7 :

 7 Figure 7: Model displacement after development of 5 rooms for base model without backfill

Figure 8 :Figure 9 :Figure 10 :

 8910 Figure 8: Model displacement after development of 5 rooms for base model with 33% backfill

Figure 9 :

 9 Figure 9: Vertical displacement (y-axis) related to backfilling degree and number of excavated rooms (x-axis)

  Backfill mechanical lab tests  Geologic and geotechnical core logs  Rock mechanical lab tests  Stress measurements and calculations  Underground support tests  Laboratory support tests for the support material as well as for the installed support elements (grouted rock bolts including the rock mass)


Christophe Auvray
email: christophe.auvray@univ-lorraine.fr

Noémie Lafrance

Mountaka Souley

Vincent Labiouse

Ageing-induced deformation of chalk rocks in underground quarries

Keywords: Chalk, ageing, underground quarries, water, strain, mechanical behaviour

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

INTRODUCTION

Alteration can be defined as a modification of the physico-chemical properties of minerals and therefore rocks, by atmospheric agents, ground waters or thermal waters [START_REF] Winkler | Applied Mineralogy, Stone properties, durability in man's[END_REF][START_REF] Foucault | Dictionnaire de géologie[END_REF]. The process is dependent on climate, water temperature, the nature of the rocks, and their degree of fissuration. Alteration generally leads to less-coherent rocks, facilitating their disintegration, rupture and eventual failure [START_REF] Kasim | An investigation of the relationship between uniaxial compressive strength and degradation for selected rock types[END_REF][START_REF] Gupta | Weathering effects on the strength and deformational behaviour of crystalline rocks under uniaxial compression state[END_REF][START_REF] Massuda | Effects of water on rock strength in a brittle regime[END_REF]. Macroscopic signs of alteration include the appearance of cracks, microfractures and macroporosity, as well as the presence of hard deposits resulting from chemical transformations [START_REF] Norbury | The description and classification of weathered rocks for engineering purposes -Geological Society Engineering Group Working Party Report[END_REF][START_REF] Chigira | Mechanism and effect of chemical weathering of sedimentary rocks[END_REF][START_REF] Oyama | Weathering rate of mudstone and tuff on old unlined tunnel walls[END_REF]. The alteration of rocks by dissolution (limestones, tuffs, sandstones, crystalline rocks), with or without mineral neoformation, has been studied by a number of authors [START_REF] Farran | L'altérabilité des roches, ses facteurs, sa prévision[END_REF][START_REF] Auger | Vieillissement par altération atmosphérique des matériaux de construction -Etude comparative in situ et en simulation[END_REF][START_REF] Furlan | Pollution atmosphérique et durabilité des pierres de construction[END_REF][START_REF] Chêne | Vieillissement accéléré de bloc de tuffeau en laboratoire sous l'effet de cycles d'imbibition-séchage[END_REF][START_REF] Chigira | Mechanism and effect of chemical weathering of sedimentary rocks[END_REF][START_REF] Gupta | Weathering effects on the strength and deformational behaviour of crystalline rocks under uniaxial compression state[END_REF].

The term 'ageing' refers to alteration of rocks as a result of human activity such as mining or quarrying. In the context of underground works, the term is used to describe the totality of mineralogical and physical modifications of the rock over time, leading to deterioration of the hydraulic or mechanical properties of the rock. The effects of time on the behaviour of an unprotected underground site consist of different phenomena responsible for degradation of the material properties: (i) degradation of constituent minerals by physico-chemical processes, leading to diminished mechanical characteristics (this is 'ageing' in its stricto sensus); (ii) deferred deformation under a constant load i.e., creep; and (iii) variations in humidity, and therefore in saturation and suction, in zones within the rock that experiences variable ventilation (these cyclical stresses can induce damage through hydromechanical coupling). The influence of each of these processes must be fully taken into account in any assessment of the temporal evolution of an underground quarry in order to provide appropriate recommendations for the solutions that should be implemented to ensure the long-term stability of a site, without endangering the safety of the public or property.

For centuries, chalk has been used as an easily-crafted and building material. The shallow underground room-and-pillar quarries where the chalk rocks were extracted have since been abandoned and now endure the effects of time and weather, potentially increasing their risk of collapse. In order to prevent this risk, INERIS, using in situ instrumentation implemented in 2004 and 2009, has been studying the behaviour of the rock in two partially flooded (subject to water table fluctuations) underground chalk quarries at Estreux and St-Martin-Le-Noeud in northern France [START_REF] Gombert | Impact du changement climatique sur la stabilité des cavités souterraines: Bilan des données de l'année[END_REF]. The two quarries, located at an average depth of 30 m, were exploited by the room and pillar method, with an extraction rate of 60-75%. The quarries were abandoned over a century ago. The original goal of the in situ characterization was to determine the impact of water-table fluctuations on their stability. These fluctuations lead to chemical interactions between the chalk and the water and affect the mechanical behaviour of the chalk due to changes in weathering conditions. Laboratory characterization was also undertaken. Nevertheless, the hydro-chemo-mechanical behaviour of the chalk remains poorly understood, despite a number of earlier studies. Due to complex interlocking phenomena acting on the mechanical behaviour of chalks, it remains necessary to address this issue in the framework of multi-scale and multiphysics approaches.

The present study is devoted to experimental characterization of the effects of ageing on the hydro-mechanical behaviour of chalk from underground quarries. The experimental tests were complemented with scanning electron microscope observations in order to identify characteristic signs of ageing in the rock.

MATERIALS AND METHODS

The Saint-Martin-le-Noeud quarry is located to the south of Beauvais, on a hillside in the town of Saint-Martin-le-Noeud. The chalk was quarried by the room-and-pillar method (overlying rock thickness ≈25m, estimated vertical stress ≈1 to 1.5 MPa), leaving galleries greater than 1.60 m in width and irregular pillars of between 2 and 4 m width. The estimated rate of extraction is 50% to 67% [START_REF] Gombert | Impact du changement climatique sur la stabilité des cavités souterraines: Bilan des données de l'année[END_REF]. The stratigraphy of Saint-Martin-le-Noeud has a near-tabular profile, with a gentle dip (≈4-10°): due to this dip, the lower parts of the quarry reach the shallow aquifer and induce several "underground lakes". The Santonian-age chalk is notable for its relatively high dolomite content (14%). The porosity -which the range of pores is tight -is on the order of 42%. The mineralogy is presented in detail in Table 1.

The Estreux quarry is situated in the town of Saint-Saulve, close to the city of Valenciennes in the North of France. Like the Saint-Martin-le-Noeud quarry, it has also been exploited by the room-and-pillar method and contains galleries of 2 to 3 m width supported by rectangular pillars measuring 1.5 to 4.5 m from side to side (rate of extraction ≈78%, overlying rock thickness ≈20 m, estimated vertical stress ≈1.8 MPa according to [START_REF] Nguyen | Influence des interactions eau-roche sur le comportement à long terme de cavités souterraines dans la craie[END_REF]. The stratigraphy has a tabular profile. The Turonian-age chalk has a relatively high glauconite content (24%), a clay mineral of the muscovite family and is not as homogenous as the dolomitic chalk in its structure.. The tight-range-of-pores porosity is around 35%. The detailed mineralogy is presented in Table 1. At each of the two sites, chalk samples were collected by drilling 120-mm-diameter, 2-to 4-mlong boreholes horizontally into the pillars. Cylindrical plugs (20 mm diameter, 40 mm height) were then taken from each of the cores at intervals of approximately 10 cm. The longitudinal axis of a plug sample is therefore parallel to the axis of the pillar (Figure 1). The samples were then used to determine the physical and mechanical properties of the chalk as a function of horizontal pillar depth. The same sampling approach has been successfully applied to gypsum in previous studies [START_REF] Auvray | The aging of gypsum in underground mines[END_REF][START_REF] Grgic | Ageing of Lorraine (France) abandoned iron mines. Rock Mechanics: a challenge for society[END_REF][START_REF] Grgic | Impact of chemical weathering on micro/macro-mechanical properties of oolithic iron ore[END_REF].

PHYSICO-MECHANICAL CHARACTERIZATIONS

Different physical and mechanical parameters were determined for all cores of up to 2 m length: rock density (dry and at saturation) ρ d / ρ sat ; porosity accessible to water, n w ; water content at saturation, W sat ; compressional ultrasonic wave velocity under saturated conditions, Vp sat , and; uniaxial compressive strength under dry / saturated conditions, RC sec / RC sat . The water used to saturate the samples was saturated (in reference to calcite) prior to the experiments in order to avoid any dissolution of the rock. The compression tests were performed using a hydraulic press of 5 kN capacity, operated at a deformation velocity corresponding to a displacement rate of 30 µm/min. Beyond the variability of results obtained from laboratory tests performed under dry and saturated conditions, we can discuss the main trends.

The results of the tests on samples from the Saint-Martin-Le-Noeud quarry are presented in Figure 2. The key findings include:

• The evolution of the four physical properties (densities, water content and porosity) clearly shows two distinct domains of variations: one defined between the wall and the depth of 75 cm, the second from the depth of 75 cm to the middle of pillar. In the first domain, wet and dry densities decreases with depth while the water-accessible porosity and water content increase. In the second domain, given the variations, a plateau can reasonably be considered for all these physical properties. • Despite showing large variation, compressional ultrasonic wave velocity generally decreases with horizontal core depth. • For mechanical characteristics (uniaxial compressive strength in saturated and dry conditions), there are also two domains of variations, with a limit of 50 cm deep against 75 cm for physical properties. In the first domain, there is no evolution of strengths: they can be considered as constant with approximately the same mean value for both states (dry and saturated). In the second domain, the dry uniaxial compressive strength increases with depth unlike the strength in saturated conditions remains on a plateau whose value is 2-3 times less than that obtained on the edge of the pillar (domain 1). Thus, for the saturated uniaxial compressive strength, the transition between the two domains is sudden. Results for the tests on samples from the Estreux quarry are presented in Figure 3. Despite systematic undulations that we do not explain:

• The evolution of all the physical and mechanical properties are monotonous except for the compressional ultrasonic wave velocity. The latter increases with depth in the first 50 cm, and decreases between 50 cm and the middle of pillar . The other physical properties decrease with distance from the outer edge of the pillar • The resistance of dry and saturated samples under uniaxial compression shows some variations, without a clear overall trend with core depth. It is nevertheless noted that the evolution of uniaxial compressive strength is identical to that shown in the second domain for Saint-Martin-Le-Noeud chalk, that is to say an increase from the wall to the middle of pillar for dry state and almost a plateau for satured uniaxial compressive strength. • Overall, the trends are clearer for the chalk of Estreux compared to that of Saint-Martin-Le-Noeud, but with lower changes.

MICROSCOPIC OBSERVATIONS

Samples from both quarries were examined by scanning electron microscopy (SEM). In samples of the Saint-Martin-Le-Noeud chalk, structural differences were observed between samples collected from the pillar-walls and those taken from the middle of the pillars. At the pillar walls, we observed euhedral multimetric-sized crystals of micro-calcite that presented no signs of neo-precipitation (Figure 4a). In contrast, at the centre of the pillars, crystals were micro-to infra-metric in size and xenomorphic (anhedral) (Figure 4e), probably due to partial dissolution. It was not possible to identify any zoning of these differences, and the decrease in crystal size and change from euhedral crystal shape both appeared to be continuous from the walls to the centres of the pillars. It was not possible to distinguish or measure the dimensions of a dissolution fringe.

In contrast to the Saint-Martin-Le-Noeud chalk, no variation in crystal size with depth was observed in the Estreux chalk pillars. However, a number of crystals in each sample plug exhibited clear dissolution marks (Figures 4g and4j). and Estreux chalks as a function of horizontal pillar depth [START_REF] Lafrance | Ageing of chalk rocks in underground quarries[END_REF].

INTERPRETATION-DISCUSSION

The effects of alteration on the physico-mechanical parameters studied here have previously been defined in samples of different rock types: charnockite [START_REF] Ramana | Quantitative studies of weathering in saprolized charnockites associated with a landslip zone at the Porthimund dam, India[END_REF], granite and marble [START_REF]Originalité des mécanismes d'altération sur les vestiges archéologiques de Délos (Cyclades, Grèce)[END_REF] and clays [START_REF] Canton | Weathering of a gypsum-calcarous mudstone under semi-arid environment at Tabernas, SE Spain: Laboratory and field-based experimental approaches[END_REF]. However, only a few studies have examined these parameters in the context of ageing (Auvray, 2003;[START_REF] Auvray | The aging of gypsum in underground mines[END_REF][START_REF] Grgic | Ageing of Lorraine (France) abandoned iron mines. Rock Mechanics: a challenge for society[END_REF][START_REF] Grgic | Impact of chemical weathering on micro/macro-mechanical properties of oolithic iron ore[END_REF]. Crystal dissolution and the observed variations in certain physical parameters should be closely related.

In the Estreux quarry chalks, the intensities of the dissolution tracks do not vary along the length of the core [START_REF] Lafrance | Ageing of chalk rocks in underground quarries[END_REF]. However, porosity accessible to water decreases from the pillar-wall to the middle of the pillar. Water (vapour and/or liquid) will progressively dissolve crystals and penetrate a pillar towards its centre. An increase in porosity accessible to water is therefore a direct effect of dissolution. This means that the amount of empty space, such as cracks and pores, should be greater at the pillar walls than in the middle of the pillars. From the SEM observations, these empty spaces appear to be in the form of dissolution cavities. The does not appear to be zoned however, and in fact, the entire pillar appears to have experienced partial dissolution. According to these hypotheses then, the resistance of the material should increase towards the centre of a pillar (the edges of the pillars being more damaged than the centres and being composed of more porous material). However, the data obtained from the uniaxial compression tests suggest that resistance is in fact highly variable, masking such a trend if it exists. It is possible that the anthromorphic changes, synonymous with ageing, are too advanced in these pillars to allow either a dissolution fringe or progressive degradation to be distinguished. Investigation of samples from a longer, horizontal borehole drilled into the rock face could shed further light on this.

In the Saint-Martin-Le-Noeud quarry samples, no dissolution marks were observed at the edges of the pillar at the scale of our observations [START_REF] Lafrance | Ageing of chalk rocks in underground quarries[END_REF]. However, marks became apparent towards the centres of the pillars. In addition, a large decrease in the size of calcite crystals was observed from the edge to the centre of the pillar. The explanation put forward for the Estreux chalk does not hold for the Saint-Martin-Le-Noeud chalk, however. Even though the density of the rock decreases from the pillar wall to its centre, the other parameters -water-accessible porosity, water content at saturation and resistance under uniaxial compression in saturated conditions -all show an overall decrease towards the centre of the pillars. In other words, the centres of the pillars are apparently more damaged than their walls. In an attempt to explain these variations, we can put forward a number of hypotheses:

• The pillar may act as a vertical drain. As a result, a greater vertical flux of water the middle of the pillar would cause massive dissolution. • Dissolution under strain would have more significant effect at the core of the pillar than at its edges. • A mechanism of aggrading neomorphism or Ostwald ripening may be involved, which would correspond to the dissolution of fine-grained crystals in the presence of an undersaturated fluid. Such dissolution would allow the simultaneous growth of larger crystals due to a phenomenon related to the surface energy of particles. Consequently, crystal sizes would be homogenized and the final size attained would be greater than it was in the initial rock. This is one of the criteria that support the aggrading neomorphism hypothesis (Morse and Casey 1988). However, with the means of observation available to us, we were unable to observe any large subhedral crystals that cover the vestiges of small rounded crystals. • A combination of the three hypotheses described above.

A dataset of measurements of resistance under compression is currently available for horizontal boreholes drilled into a phosphate chalk quarry (Malogne quarry at Cuesmes), near Mons in Belgium [START_REF] Funcken | Procès verbal-Stabilité de la carrière de la Malogne[END_REF]. The resistance values measured in the Malogne chalk appear to decrease towards the interior of the pillar (Figure 5), thus supporting the findings of our study and allowing us to favour a hypothesis of preferential vertical drainage down the middle of a pillar. We can therefore infer that ageing of the chalk also occurs at Estreux, but it is not possible to identify the mechanism or mechanisms responsible.

CONCLUSIONS

At the quarry at Saint-Martin-Le-Noeud, we have observed evidence for temporal changes to the chalk due to human activity, changes which can therefore be defined as ageing. The intensity of dissolution marks in calcite grains appears to decrease outwards, from the inner-core of the pillar towards the pillar walls. However, rock density decreases in the opposite direction, from the wall towards the centre, and parameters such as water-accessible porosity, water content at saturation and resistance to uniaxial compression under saturated conditions also present an overall decrease from the wall towards the middle of the pillars. Four hypotheses can be put forward to explain the observed variations: (i) vertical drainage that is greater in the interior of the pillar than at its edges; (ii) a higher degree of dissolution under strain down the middle of the pillar; (iii) aggrading neomorphism or Ostwald ripening; (iv) a combination of the three.

At the site at Estreux, ageing of the glauconitic chalk already appears to be at an advanced stage. This may explain why we were unable to distinguish zoning of either dissolution marks or resistance to uniaxial compression. To test this hypothesis, a long horizontal borehole could be drilled in the chalk face, thus allowing us to investigate the nature of the chalk at depths beyond those regions affected by ageing and to reach areas of intact 'healthy' chalk.

The paralleling of SEM analysis and various physico-mechanical properties has permitted us to probe into the ageing-induced changes of different chalks with efficiency and little material. However, a certain amount of time is required just for preparing the samples which might increase because of rocks variability. It may not be wise to apply this method in too much non homogenous rocks, especially if a global geological study (dip, thickness of the geological disposal, global rheology…) was not made prior. From a broader perspective, the distribution of mechanical stresses inside the pillar must be raised in parallel of this method, especially when dealing with soft, complex and as easily deformable as chalk material. It may therefore provide a part of explanation to the variations encountered in the decrease of the parameters from wall to centre of the pillar.

Figure 1 -

 1 Figure 1 -Sampling method.

Figure 2 -Figure 3 -

 23 Figure 2 -Variations in physico-mechanical parameters with horizontal core-depth, Saint-Martin-Le-Noeud quarry.

Figure 4 -

 4 Figure 4 -Structural and morphological variations ( : dissolution marks) in the Saint-Martin-Le-Noeud and Estreux chalks as a function of horizontal pillar depth (Lafrance et al., 2014).

Figure 5 -

 5 Figure5-Variations in uniaxial compressive strength (natural water content) as a function of horizontal pillar depth[START_REF] Funcken | Procès verbal-Stabilité de la carrière de la Malogne[END_REF] Mons quarry

Table 1 -

 1 Mineralogical compositions of chalk from the two quarries. (Spa: sparite; Mic: micrite; Qtz: quartz; Glau: glauconite; Dol: dolomite; Opq: opaque minerals)

		Spa	Mic	Qtz	Glau	Dol	Opq
		(%)	(%)	(%)	(%)	(%)	(%)
	Saint-Martin	17	68	0	0	14	1
	Estreux	24	50	2	23	0	1

ACKNOWLEDGMENTS

We would like to thank Mr Pierret (Infrastructure Department at INFRABEL) and Mr Funcken (Direction Générale des routes et Bâtiments du Service Public de Wallonie) for their assistance in this study. We also want to thank the Conservatoire Naturel d'Espaces de Picardie for authorization to access to the quarry, major site wintering chiroptera.


Mountaka Souley
email: mountaka.souley@ineris.fr

Mehdi Ghoreychi

Jean-Bernard Kazmierczak

Gilles Armand

D Seyedi

J.-B Kazmierczak

An anisotropic elastoplastic model for geomaterials and numerical implementation

Keywords: Constitutive model, structural anisotropy, induced anisotropy, numerical implementation, verification

INTRODUCTION

Numerous experimental results available in the literature indicate that most sedimentary and metamorphic rocks, such as shales and slates, display a strong anisotropy of behavior and strength. These types of rocks usually exhibit some preferentially oriented structures or possess distinct bedding planes, which results in transversely isotropic behavior at the macroscopic scale: this is particularly the case of Tournemire shale [START_REF] Niandou | Laboratory investigation of the mechanical behavior of tournemire shale[END_REF] or the Callovo-Oxfordian claystone [START_REF] David | Anisotropy of elastic, magnetic and microstructural properties of the Callovo-Oxfordian shales[END_REF].

Based on the experimental results, various failure criteria for anisotropic materials have been proposed where a large review can be found in [START_REF] Duveau | Assessment of some failure criteria for strongly anisotropic materials[END_REF]. Furthermore, various constitutive models were developed in order to approach the anisotropic mechanical behavior. These are mainly (a) the empirical models based on the theory of variational cohesion and / or friction (e.g., [START_REF] Mclamore | The mechanical behavior of anisotropic sedimentary rocks[END_REF][START_REF] Saroglou | A modified Hoek-Brown failure criterion for anisotropic intact rock[END_REF][START_REF] Wang | Three-dimensional shakedown solutions for anisotropic cohesive frictional materials under moving surface loads[END_REF]; (b) the models built on the concept of ubiquitous joint with several planes of weakness, (c) the models where damage and/or plasticity are incorporated and formulated in the framework of irreversible thermodynamics [START_REF] Pietruszczak | Modeling of inherent anisotropy in sedimentary rocks[END_REF][START_REF] Yu | A Transversely Isotropic Damage Model for Boom Clay[END_REF].

Since 2000, the French National Radioactive Waste Management Agency (Andra) has been constructing an Underground Research Laboratory (URL) at Bure with intent to demonstrate the feasibility of a geological repository in the Callovo-Oxfordian claystone (COx) formation by collecting in situ experimental data. The excavation of galleries at the main level of the laboratory showed a significant fracturing induced by the excavation [START_REF] Armand | Geometry and properties of the excavation induced fractures at the Meuse/Haute-Marne URL drifts[END_REF] in addition to structural or inherent anisotropy observed during the laboratory tests.

Alongside this in situ investigation program and the laboratory characterization, several developments of rheological models and their implementation in numerical codes have also been undertaken. A macroscopic isotropic visco-elastoplastic model, which accounts for the impact of induced damage and fracturing on the delayed strain rates, has already been proposed in our previous studies [START_REF] Souley | Modelling of the viscoplastic behaviour including damage for deep argillaceous rocks[END_REF]. However, this model did not allow a full explanation of the different types of observed fractures, their extension, some convergence anisotropy and ratios, or the reproduction of some in situ pore pressure observations. For this purpose a macroscopic anisotropic phenomenological model (including both the structural anisotropy and the induced anisotropic plasticity) has been proposed in this paper. It consists in an extension of the implicit model for strength anisotropy implemented in ITASCA codes. The basic assumption is that the failure of an anisotropic material is due to either fracturing of bedding planes (or weakness planes) or the failure of the rock matrix. The formalism is therefore based on approaches with weakness planes commonly called "Discontinuous weakness plane" or "Ubquituous joints" (aimed to account for rock mass strength and anisotropy within large-scale continuum models, [START_REF] Sainsbury | Simulation of rock mass strength anisotropy and scale effects using a Ubiquitous Joint Rock Mass (UJRM) model[END_REF] with one of the first models was proposed by [START_REF] Jaeger | Shear failure of anisotropic rocks[END_REF] and known as "the single plane of weakness theory".

BRIEF REVIEW OF THE ANISOTROPIC BEHAVIOR OF GEOMATERIALS

In the literature, many experimental studies have been conducted to study the mechanical behavior of anisotropic geomaterials characterized by a large directional dependency: existence of an anisotropy of the elastic modules and the strengths as a function of the load orientation. These results were also confirmed by the study of several materials having generally a transverse isotropy by Attwell & Sandford (1974), [START_REF] Niandou | Laboratory investigation of the mechanical behavior of tournemire shale[END_REF], with the main features: (a) the strength of anisotropic material depends on the stress state and the loading orientation with respect to the stratification, (b) the minimum strength is generally reached when the angle of inclination β between the major stress and the weak planes ranges from 30 to 60°, (c) the ratio of resistance in the two main directions related to the anisotropy or the ratio between the minor and major strengths are used to quantify the degree of anisotropy, (d) failure may occur in very complex ways, combining sliding and separation along the planes of weakness, and a shear band in the matrix.

Research on Callovo Oxfordian claystone by measurement of the propagation velocities of ultrasonic waves showed a structural anisotropy. The work conducted by [START_REF] David | Anisotropy of elastic, magnetic and microstructural properties of the Callovo-Oxfordian shales[END_REF] confirms the anisotropy of the wave propagation speed. The anisotropy of both the wave propagation speed and the dynamic modules has been also observed on the mechanical laboratory tests with acoustic emissions monitoring [START_REF] Sarout | Anisotropy of elastic wave velocities in deformed shales. Part I: Experimental results[END_REF].

The anisotropy of the mechanical behavior of the Callovo Oxfordian was studied by a series of tests carried out in different orientations with respect to the stratification. The results show an anisotropy of deformation modulus, with a modulus ratio II E E ⊥ between 1.05 and 1.4.

Traditional measures of modules (from compression tests) show that the anisotropy could be larger, but still less than a ratio of 2. The ratio between the strength as a function of β and the strength in the isotropic plane varies between 1 and 0.8. Results showed that the strengths in the directions parallel and perpendicular to the stratification are close and slightly higher than those measured at 45 ° with respect to the bedding plane with a ratio close to 0.8 for the majority of the performed tests. However, these results need to be confirmed by additional measures according to different mineralogy.

Models based on the concept of "Discontinuous Weakness plane" allow a description of the physical mechanism related to the process of failure. In general, an anisotropic material is viewed as an isotropic body containing some planes of weakness at the macroscopic scale. In this case the basic assumption is that the failure of the anisotropic material is due to the fracture planes as well as to the rock matrix. Assuming that the weakness planes are oriented along the cracks and based on the extension of Griffith theory, several criteria are proposed in the literature [START_REF] Walsh | A fracture criterion for brittle anisotropic rock[END_REF][START_REF] Hoek | Empirical strength criterion for rock masses[END_REF][START_REF] Hoek | Strength of jointed rock masses[END_REF]. Generally speaking, the rock matrix is modeled by nonlinear criteria involving the three invariants of stresses, while sliding along the planes of weakness follows the Mohr-Coulomb criterion or the JRC-JCS model of [START_REF] Bandis | Application of a New Numerical Model of Joint Behaviour to Rock Mechanics Problems[END_REF] developed for rock joints [START_REF] Wang | A constitutive model for the deformation of a rock mass containing sets of ubiquitous joints[END_REF].

The proposed model is based on the concept of ubiquitous joint and yielding may occur either in the solid or along the weak planes depending on the stress state, the orientation of the weak planes and the material properties of solid and weakness planes.

FORMULATION OF THE MODEL

In the absence of sufficient data (especially quantitative) on directional dependency of the strengths and plastic flows (hardening/softening), we firstly propose a relatively simple model that can integrate the anisotropy of the elastic behavior and that of rupture, but sufficiently generic and phenomenological so that it can be subsequently enriched on the basis of (a) new knowledge and experimental data, (b) other phenomena not taken into account here (damage in the framework of continuum damage mechanics, hydromechanical couplings, time-dependent behavior, etc.).

Case of the intact rock matrix

The rheological model proposed is based on results obtained from different laboratory characterizations. The main features of the short-term mechanical behavior observed on the samples of claystone under triaxial tests can be summarized as follows: the behavior is linear under low deviatoric stresses; the loss of linearity of stress-lateral strain curves begins approximately at 50% of the peak value of the deviator stress. Under low confining pressures, a brittle failure of the samples is observed and corresponds to the formation of a shear band inclined with respect to the sample axis. There is a strong dependence of the mechanical behavior on the confining pressure, marked by a transition from a brittle toward a ductile behavior.

It has been shown that a failure criterion based on the Hoek-Brown ( 1980) criterion (eq. 1) is suitable for describing the shear strength of the COx claystone and as a result an elastoviscoplastic and isotropic model was proposed [START_REF] Souley | Modelling of the viscoplastic behaviour including damage for deep argillaceous rocks[END_REF]. However, since the COx claystones in their natural state have no tectonic fractures, classical signification of the Hoek & Brown parameters cannot be used for determining the rock mass strength from the results of tests on samples (and taking into account the fracturing of rock).

2 3 3 c 3 c m s σ = σ --σ σ + σ (1)
where σ 1 and σ 3 are major and minor (absolute value) principal stresses (compressive stresses are negative), m and s are material constants, σ c is the value of stress reached at the initiation and peak under unconfined compressive condition.

After some manipulations, the generalization of equation ( 1) in the space of three stress invariants (p, q, θ) leads to the following criterion for damage initiation and peak strengths:

2 m 2 s c 4 cos cos sin B F q t q p 3 A 3 A 3   θ θ θ = + - + -     (2)
where p represents the mean stress, q the generalized deviatoric stress and θ the Lode's angle, A and B are two independent parameters (A = m σ c and B = s σ c 2 , evaluated at the initiation and peak). The parameter t c (≤ 1) has been introduced in order to account for a potential difference of material strength for compression and extension loading paths.

Under tensile stresses, we propose a «cutoff» of mean stress to the value of tensile strength σ t (in triaxial condition, which is by default: B/A), this leads to following yield surface in tension:

m t t F p = -σ (3)
As shown by the experimental results, we assume that the behavior of COx claystone is transversely isotropic. Assuming that the results (small changes in strength) of tests previously discussed or Brazilian tests are representative of the strength characteristics of the rock matrix, then it follows that a quasi-isotropic failure can be considered as a first approximation of rock matrix failure. The induced anisotropy due to the change of stresses will be approximated by an assumption of the existence of a weakness plane with a certain orientation. As a result, the independent parameters A and B do not change with the loading direction.

In 2005, some complementary tests carried out on claystone samples made it possible to specify the evolution of the irreversible volumetric deformation. From these tests, it came out that the hypothesis of associated flow rule tended to over-estimate the material dilatancy. Indeed, experimental data [START_REF] Chiarelli | Modelling of elastic-plastic damage behaviour of a claystone[END_REF][START_REF] Zhang | Experimental study of the hydro-mechanical behavior of the Callovo-Oxfordian argillite[END_REF]) exhibit generally a contracting behavior. At best, a dilatancy develops in the last phases of the tests according to the confining pressure. The following form of plastic potential is adopted.

( )

m s c cos sin G t q p 3 3   θ θ = - + β γ     (4)
The rate of dilatancy β(γ), where γ (the internal flow variable) is the plastic distortion, varies between a minimum value β 0 and a maximum value β m , its evolution was first proposed by [START_REF] Chiarelli | Modelling of elastic-plastic damage behaviour of a claystone[END_REF] and used by many authors thereafter:

( ) ( ) b m m 0 e β -γ β γ = β -β -β (5)
where b β is the plastic flow velocity.

Finally an associated flow rule is assumed in tension.

For residual strength, the same shapes of yield function as the peak are suitable when the confining pressure does not exceed the transition stress between the ductile and softening behaviors, σ 3 bd . Beyond this confining pressure, residual strength coincides with the peak yield function. The needed parameters of residual strength (when σ 3 <σ 3 bd ) are s r and m r . By default s r is taken equal to 0 and m r is back-calculated so that the residual strength intercepts the maximum resistance for σ 3 bd of confining pressure.

The strain hardening in pre-peak and strain softening in post peak are modeled as a linear change of A with respect to the internal plastic variable (from the initiation A i to the peak A p ). We also considered a parabolic evolution of parameter B as a function of γ, with horizontal tangents at peak and the beginning of residual phase, respectively, for the hardening in pre-peak and softening in post-peak.

Case of weakness planes (induced anisotropy)

The failure criterion on the weak plane, whose orientation will be defined later, is a composite Mohr-Coulomb envelope with strength «cutoff» in tension (similar to the FLAC 3D Ubiquitous Joint Model). This is justified since we can associate each plane of weakness as a dilatant (or not) rock joint (depending on the roughness type: primary and / or secondary). At this stage, we simply consider a basic failure criterion. The projection of a state of stress on the composite envelope is controlled by a non associated flow and an associated flow rules, respectively for shear and tensile failure.

Let (s, t, n) be the local system of reference axes related to the weakness plane consisting of (s, t) and n pointing in the direction of the unit normal of this plane. The magnitude of the tangential (shear) stress component on the weak plane, and the associated shear strain variable are:

The generalized stress vector used to describe weak-plane failure has four components: σ ss , σ tt , σ nn , τ. The components of the corresponding generalized strain vector are: ε ss , ε tt , ε nn , γ wp . The weak-plane failure criterion with a tension cutoff is expressed in terms of (σ nn , τ) as follow:

wp s nn wp wp wp t t nn wp F tan C F  = τ + σ φ -   = σ -σ   (7)
where φ wp , C wp and t wp σ are respectively the friction, cohesion and tensile strength of the weak plane.

The potential function is composed of two functions used to define shear (non associated) and tensile (associated) plastic flow, respectively:

wp s nn wp wp t nn G tan G  = τ + σ ψ   = σ   (8)
where ψ wp is the weak-plane dilation angle.

Several experimental studies have established a correlation between the mechanical behavior of geomaterials and their microstructure [START_REF] Kranz | Microcracks in rocks -a review[END_REF]. Microscopic observations on samples of several rocks subjected to different stress levels showed that microcracks follow the deformations of extension under certain preferential directions. This phenomenon induces a dilatancy and induced anisotropy. Based on different experimental results, some authors propose that the driving force of the induced anisotropic damage could be directly related to the positive part (extension) of the total strain tensor after the spectral decomposition introduced by Orizt (1985) and supplemented by [START_REF] Ju | On energy-based coupled elastoplastic damage theories: constitutive modelling and computational aspects[END_REF].

As the failure along the plane of weakness may be either in shear or in tension, it can be supposed that the weakness plane remains normal to the direction of the minor principal stress. Indeed, the fact to consider that the plane of weakness is related to the positive part of the strain tends to exclude the shear failure that may occur under a compressive loading path. In the numerical implementation, the orientation of the weakness plane (in each element) is arbitrary and corresponds to an input data to be introduced by the user. By default, the normal of weakness plane is taken as the direction of the minor principal stress.

Constitutive equations and numerical implementation

Assuming that only small strains occurred, the total strain increment dε can be subdivided in an elastic part e dε and a plastic part p dε :

e p d d d ε = ε + ε with p G d ∂ ε = λ ∂σ (9)
where λ is the plastic multiplier and G is the plastic potential for shear or tensile failure, for both solid rock and weakness plane related to the induced anisotropy.

If C represents the fourth order elastic compliance tensor for transversely isotropic geomaterial, the incremental expression of Hooke's law in terms of generalized stress and strain tensors has the form:

e d Cd σ = ε (10)
Finally, the consistency condition dF(κ, σ) = 0 allows to express the stress increment as a function of total strain increment:

F G C : C : d C d F G F G : C :       ∂ ∂ ⊗       ∂σ ∂σ       σ = - ε   ∂ ∂ ∂ ∂ -   ∂σ ∂σ ∂κ ∂ϖ     (11)
where F is the yield function; κ is the accumulated plastic strain (e.g. γ for matrix failure or γ wp for failure along the weakness plane); ϖ is the generalized deviatoric stress (q) for matrix failure and the tangential stress (τ) for failure along the weakness plane.

For both (rock matrix and weakness planes) failures, a function representing the diagonal between F s = 0 and F t = 0 is evaluated in order to select the type of failure (in shear or in tension). At each step, the computation of new stresses is achieved by testing firstly failure in the solid matrix, then along the plane of weakness. More precisely, the main procedure is summarized below.

• The first approximation of stress tensor σ i , is evaluated by adding to the previous stress tensor the stress increment computed from the total strain increments and the Hooke's law. • Computation of the yield function for rock matrix, F m (p i , q i , θ i ) according to eqns 2 and 3. If the stress state, σ i , satisfies the yield criterion: F m (p i , q i , θ i ) > 0, the new increment of stresses is computed as well as the new stress tensor, σ o . If F m (p i , q i , θ i ) ≤ 0, the current stress components are:

σ o = σ i .
The resulting stress tensor, σ o , is then examined for failure on the weak plane.

• The corresponding stress components in the local axes, σ' o , are computed using the transformation from global x-, y-, z-axes to the local system of reference s-, t-, n-axes. • Computation of the yield function for weak plane, F wp (σ nn ' o , τ' o ) according to equation ( 6). If the stress state σ' o verifies the yield criterion F wp (σ nn ' o , τ' o ) > 0, the new increment of stresses is computed as well as the new stress tensor, σ n . If F wp (σ nn ' o , τ' o ) ≤ 0, the current new stresses are:

σ n = σ' o .
• Back to global system of reference x-, y-, z-axes.

• Parameters depending on the internal plastic variable are updated. This routine has been written in C++ and compiled as a Dynamic Link Library (DLL) file that can be loaded whenever it is needed.

All parameters have already been identified from triaxial tests for the short-term response [START_REF] Su | Constitutive Models for the Meuse/Haute-Marne Argillites -MODEX-REP, European Commission -Nuclear science and technology[END_REF][START_REF] Chiarelli | Modelling of elastic-plastic damage behaviour of a claystone[END_REF]. Note that, results of laboratory tests showed a wide variability of peak and residual plastic strains, they can be viewed as calibration parameters for in situ applications. Parameters related to the failure along planes of weakness were not for the time based on a characterization campaign: the choice of the considered values is arbitrary and deserves further attention in future studies. Also for reasons of simplicity the elastic characteristics along transverse planes of weakness are considered identical to those of the matrix as by default (E 1 =E s , Young modulus in the isotropic (matrix) and weakness plane; E 3 =E n , Young modulus in the direction normal to the isotropic (matrix) and weakness plane; G 13 =G sn , shear modulus in the anisotropic plane; ν 12 =ν st and ν 13 =ν sn , Poisson coefficient respectively in the isotropic and anisotropic planes). Table 1

-0.1 β m 0.3 b β 600 t c 1 (a)
G 13 =1807 MPa is computer based on the relation of Lekhnitski (1981)

VERIFICATION AND APPLICATION

Triaxial compression path

In order to verify the model (constitutive equations and numerical implementation), several triaxial compression tests with confining pressures of 0.01, 1, 5, 10, 20, and 30 MPa are simulated. Figure 1 presents the deviatoric stress -axial, lateral and volumetric strain curves for different confining pressures. The difference between predicted strengths and those obtained in laboratory tests is related to the lower values of the parameters s i , σ ci , m i , s p , σ cp , and m p used in this study. The lower values are appropriated for modeling of in situ structures, whereas the mean values allow to well reproduce the laboratory strengths. From figure 1a, we note that the post-peak behavior is confining pressure dependent: the transition stress between brittle failure and ductile behavior is clearly marked and with an input value of transition stress, σ 3 bd , about 25 MPa. The resulting curves display four regions (elastic, hardening in pre-peak, softening in post-peak and residual phases) when the confining pressure is below the transition stress level σ 3 bd . It displays three regions (elastic, hardening and perfect plastic phase) under high confining pressure. These curves are qualitatively similar to the experimental ones (which were not reported herein).

Figure 2 shows a comparison in terms of the elastic limit, the peak and residual strengths between the simulation (corresponding values in Figure 1) and theory (Eq. 1). A good agreement can be observed as numerical and analytical solutions coincide. More precisely, the relative error for peak strength is less than 0.1%, and 3% for the elastic limit and residual strength respectively. This validates the numerical implementation of the proposed anisotropic elastoplastic model in FLAC 3D . The aim of this section is to provide a verification of the implementation for non-triaxial stress paths and to show numerically the ability of the implemented model to evaluate the extent of damaged (in pre-peak) and failed (post-peak) zones around a circular underground excavation. The GCS drift excavated in the direction of the principal horizontal major stress is selected for this application. This is one of the several mine-by experiments set up at the main level -490 m of the Meuse/Haute-Marne URL where an anisotropic in situ stress state has been measured [START_REF] Wileveau | Complete in situ stress determination in an argillite sedimentary formation[END_REF]. The GCS drift has a circular section with a 2.6 m radius and can be considered as an isolated gallery.

We deliberately chose a drift oriented in the direction of the principal horizontal major stress, leaving a very low stress anisotropy ratio (between the principal minor horizontal stress and principal vertical stress) of 0.98 in the drift section and for which earlier simulations performed in the framework of the isotropic continuous approaches failed to reproduce the extension of fractured zones and the ratio of convergence observed in situ.

Taking into account the knowledge on the in situ stress state, σ v , σ h and σ H are considered to be equal respectively to 12.7, 12.4 and 16.12 MPa, corresponding to the maximal stress ratio of 1.3. In the numerical model centred at the depth of 490 m, the three initial in situ stresses are distributed according to the depth. Roller displacement conditions were applied to the base of the model whereas boundary stresses (depending on the depth) were applied to the others planes (top and lateral limits of the model).

Figure 3 shows the extent of plastic zones around the circular opening. In this figure the legend "hard", "soft", "res" correspond respectively to the shear failure of the rock matrix in the pre, post peak and residual region. "Shear" and "tens" are associated to the shear and tensile failures along the plane of weakness. Plastics radius in the vertical and horizontal directions are respectively 3.4 and 5.3 m, that is to say 0.65 D and 1.02 D. These extensions are in the range of fractured (connected or not) zones extensions observed around the drifts oriented in the direction of the major horizontal stress [START_REF] Armand | Geometry and properties of the excavation induced fractures at the Meuse/Haute-Marne URL drifts[END_REF].

For this instantaneous mechanical calculation, the horizontal and vertical convergences are respectively 5.2 and 4 cm, this corresponds to an initial convergence ratio of 1.3. Figure 4 shows the profiles of radial, orthoradial and axial stresses along two radial lines at approximately 0° and 90° with respect to x-axis as a function of the radial distance from the gallery center. From the profile of orthoradial and axial stresses, one can distinguish different regions (elastic, pre-peak and postpeak and residual, without particular transition between the post-peak and the residual domains) through the slopes of curves. In particular at θ = 90°, the first loss of slope can be noticed along the profile of orthoradial and axial stresses: the corresponding radial distance (approximately 3.4 m) is in accordance with the previous investigation of plastic zone extent. Between this radial distance and 6 m, there is a nonmonotonic evolution of orthoradial and axial stresses in relation to the plastic area formed at ± 65 ° with respect to horizontal direction.

CONCLUSION

This paper presents the development of an anisotropic constitutive model for transversely isotropic geomaterials, its numerical implementation in the three-dimensional code FLAC 3D , as well as its verification and validation. Firstly, a non-linear elastoplastic model based on the interpretation of laboratory test performed using the Hoek-Brown failure criterion, is presented. The elastic anisotropy, the almost isotropy of strength observed on the available tests carried out in several directions, the induced anisotropy observed on the in situ structures are the main ingredients for this model. The rock behavior is characterized by a non-linear behavior in pre peak, a non-linear softening in post-peak and a perfectly plastic behavior in the residual phase. As observed on laboratory tests, the model takes into account the transition between the softening in post-peak under low confining pressures and the ductile behavior under high confining pressures.

Simulation of triaxial compression tests at different levels of confining pressure provides a verification of the implemented model. The resulting curves display four regions (elastic, hardening in prepeak, softening in post-peak and residual phase) when the confining pressure is below the transition stress, and three regions (elastic, damage and perfect plastic phase) under high confining pressure. In addition, the elastic limit, the peak and residual strengths derived from these simulations are compared with the theoretical envelopes: the corresponding relative error does not exceed 1%.

The ability of the implemented model to reproduce the plastic zones around an underground drift excavated in the Callovo-Oxfordian claystone is successfully tested. The results provide new insights on the understanding of the deformation mechanisms observed around the structures of Meuse/Haute-Marne Underground Research Laboratory.

Next step will be to enrich the proposed model on the basis of new knowledge and experimental data, and/or other phenomena that have not been taken into account here (damage in the framework of damage mechanics, hydromechanical couplings, time-dependent behavior etc.).

Figure 1 .

 1 Figure 1. Triaxial compression tests: (a) numerical simulations, stress vs strains; (b) laboratory tests, stress vs axial displacement

Figure 2 .

 2 Figure 2. Elastic limit, peak and residual strengths: numerical and analytical solutions Application to a circular opening of Meuse-Haute Marne Underground Research Laboratory

Figure 3 .

 3 Figure 3. Extent of plastic zones around the opening

Figure 4 .

 4 Figure 4. Radial, orthoradial, axial and shear stresses along radial lines (a) at 0° (b) at 90°

  provides the list of parameters and input values.

				Table 1. Values of input parameters
	Parameter Value	Parameter Value	Parameter Value	Parameter Value
	E 1 =E s	5600 MPa m i	1.5	m p	2	σ 3 bd	25 MPa
	E 3 =E n G 13 =G sn ν 12 =ν st ν 13 =ν sn	4000 MPa s i (a) σ ci 0.3 C wp 0.25 σ wp t	1 9.6 MPa 1 MPa 1 MPa	s p σ cp φ wp ψ wp	0.128 33.5 MPa 15° 5°	0 m r (output) 2.172 s r γ peak 3.75 10 -3 γ residual 8.5 10 -3
	β 0


Emanuela Matrullo

Isabelle Contrucci

INERIS), I. Contrucci (INERIS Pascal Dominique

Mustapha Bennani

Hideo Aochi

Jannes Kinscher

Pascal Bernard

Pascal Bigarre

Analysis and interpretation of induced micro-seismicity by flooding of the Gardanne coal basin (Provence -Southern France)

and research institutions in France or abroad, or from public or private research centers.

Analysis and Interpretation of Induced Microseismicity by Flooding of the Gardanne Coal Basin (Provence -Southern France)

Introduction

Mining activity in Provence has left many underground voids that can cause land subsidence or collapses in areas of high population density. Part of these underground voids are now flooded and causes seismicity, whose strongest events were felt by the population. To mitigate these risks, a realtime microseismic monitoring was implemented in sensitive areas. And since the cessation of mining concessions in Provence, INERIS (Institut national de l'environnement industriel et des risques) has installed 5 permanent microseismic monitoring networks in the Gardanne coalfield which was operated by Charbonnages de France (Figure 1). Then, BRGM-DPSM (Bureau de recherches géologiques et minières -Département Prévention et sécurité minère) on behalf of the State, asked INERIS to continue the microseismic monitoring.

The flooding of mine voids and fluctuations in the groundwater level induced many seismic events that have been recorded since the beginning of operational monitoring. An example is the seismic crisis of November 2012, where magnitude events > 2.5 were felt by the population [START_REF] Dominique | Analyse préliminaire de la crise sismique de Fuveau (Bouches-du-Rhones)[END_REF][START_REF] Contrucci | Activité microsismique et caractérisation de la détectabilité des réseaux de surveillance du bassin houiller de Gardanne[END_REF][START_REF] Bennani | Note d'expertise complémentaire relative aux événements microsismiques ressentis en surface dans la zone de Fuveau -Gréasque, dans le bassin minier de Gardanne[END_REF]. In 2013, 5 temporary seismic stations were deployed by BRGM-DPSM specifically to monitor the activity of that particular area. Early 2014, due to an exceptionally high rainfall and achieving maximum pumping capacity in the mining area, the groundwater level of the mining reservoir rises rapidly. To anticipate the possibility of a renewed seismic activity, a temporary seismic network of 5 stations was deployed by INERIS, to complete the coverage of the existing network. Thus, the coal basin Gardanne is now monitored by five permanent microseismic stations and 10 temporary seismic stations (Figure 1). In December 2014, a new microseismic crisis was felt by people leaving in the area.

Scientific objectives

The flooding of the mine generates significant seismicity highlighted upon cessation of water pumping and observed since the beginning of the operational monitoring of the basin, confirming a high reactivity of mining voids. Strong seismic events (> 2.2) were recorded raising serious concerns of the people (Figure 1a). The question of events induced by the mine and flooding, or triggered events with fault re-activation, or natural events is asked. Depending on the response, the estimate of the likely maximum magnitude in the mining area may differ significantly. In this context, the objective is to engage in-depth analyses by taking advantage of the temporary network to understand the mechanisms of the observed seismicity. This will include clarifying the relationships between seismicity and rising waters in old mine working, geology and local tectonics.

We will focus on the spatio-temporal distribution of the seismic events which allows to monitor the dynamics of faulting and to estimate the extension of the fracture system. The determination of the mode of rupture (focal mechanisms) of the events is also a very important information to characterize the stability of the abandoned mine. The present work aims at a better understanding of the mechanical characteristics of the fracture system based on the observation, in particular on an accurate analysis of the low magnitude seismicity. The general aim of the work will be pursued, according to the availability of the dataset from high quality recording of the dense arrays. We will focus on the characterization of the microseismicity in terms of the accurate hypocentral locations, local and moment magnitudes and accurate focal mechanisms determination.

Methods

Such understanding of the faulting or underground mine collapse mechanisms needs accurate seismic event locations and precise determination of the focal mechanisms. In order to obtain reliable event locations and focal mechanisms of the low magnitude seismicity, an accurate 1D velocity model of the structure at depth is determined solving the coupled hypocenter-velocity inverse problem for more than 200 local earthquakes recorded. The station corrections reflect the inhomogeneous near-surface features.

In addition we tested and validated different methods of source localization, based on various data parameters (travel time, polarization, and combinations) using probabilistic approach [START_REF] Contrucci | Management of post-mining large-scale ground failures : blast swarms field experiment for calibration of permanent microseismic early-warning systems[END_REF][START_REF] Lomax | Probabilistic earthquake location in 3D and layered models: Introduction of a Metropolis-Gibbs method and comparison with linear locations[END_REF].

Spectral amplitude ratios will be analyzed to provide information on the types of focal mechanisms (tensile versus double-couple, Kinscher, 2015). The application of a general automatic moment tensor inversion approach [START_REF] Dahm | Automated moment tensor inversion to estimate source mechanisms of hydraulically induced micro-seismicity in salt rock[END_REF], Sen et al., 2013, e.g.) is very difficult for the high frequency microseismic data record because of the great complexity of the seismic signals. The source mechanism will be constrained by using an inversion approach considering the observed amplitude spectra of different stations in order to explore the best source model with respect to one single event (Kinscher, 2015). Spectral amplitude will be also used to determine the source seismic energies (seismic moment and magnitude) in order to calibrate a robust magnitude law. The data processing will be automated and applied at the whole dataset of events in order to obtain a complete catalog describing the evolution of the mine seismicity. Particular attention will be paid on the two principal crises happened in November 2012 and December 2014 (Figure 1a). These two micro-earthquake swarming sequences repeatedly occur in specific sub-areas and appear influenced by the fluid migration. For this purpose, the distribution and space-time migration of the micro-seismicity will be analyzed in terms of variations in seismic wave propagation properties. In particular to better understand the role of the fluids in the fault failure process the ratio between compressional-wave and shear-wave velocity (Vp/Vs) will be analyzed.

Conclusions

Some preliminary analysis on the Vp/Vs ratio and P-wave velocity model was performed. We found a high Vp/Vs (1.85) and a 1D velocity model consistent with the lithologies present in the investigated area. Very small station residuals (between -0.05 and 0.05) are apparently related to the geometry of the mine. The analysis reveals and confirms the presence in the studied area of a carbonated layer partially fractured and permeated of water. The seismicity appears correlated with the waterfront. This suggest the high influence of the water migration on the micro-earthquakes triggering (Figure 2). The analysis on the source characterization will give more information about the fracture process.

Figure 1

 1 Figure 1 (a) Map of the recent seismicity (2008-2014) detected and located by INERIS. The colours of the circle are related to the duration magnitude taken from the catalogue. Several symbols are used for the different networks present in the area. The yellow circle represents the area of high population density of 400m around the stations. The limit of the waterfront is also represented with a blue line. The red zone represents the extension of the "Grande mine" mine working layer. (b) INERIS permanent network configuration.

Figure 2

 2 Figure 2 Correlation of the number of events detected by the BRGM temporary network on the crisis area (see red star on Fig. 1) (top figure) and (1) the daily rainfall (green), (2) the water level (red) and (3) the velocity of the water level change (blue) for the period July 2013-August 2014 (bottom figure). The daily rainfall is provided by METEO FRANCE, the number of events and the water level data are provided by BRGM-DPSM.

(a) (b)