Claudio Lobos

Marek Bucki

Nancy Hitschfeld

Yohan Payan

Mixed-element mesh for an intra-operative modeling of the brain tumor extraction

Keywords: Modified-Octree, Mixed-elements, Finite Elements, Real-time Application, Patient-specific Simulation, Region Of Interest

This paper presents a modified-octree technique that generates a mixedelement mesh. The final output mesh consider cubes, prisms, pyramids and tetrahedra. This technique is optimized for brain tumor extraction simulation in a real-time application. The proposed technique is based on the octree algorithm with a specific constraint: elements will be split only if they intersects a certain region of interest. With this approach we pursued a refined mesh only in the path from the skull opening point to the tumor. Fast computation by the Finite Element Method (FEM) is achieve thanks to the local refinement. Examples are given and comparison with other approaches are presented.

Clinical Background

Accurate localization of the target is essential to reduce morbidity during a brain tumor removal intervention. Image guided neurosurgery nowadays faces an important issue for large skull openings, with intra-operative changes that remain largely unsolved.

Once the skull is open a deformation of the brain naturally occurs. This phenomena is known as "the brain-shift". The causes of this deformation can be grouped by:

• physical changes (dura opening, gravity, loss of cerebrospinal fluid, actions of the neurosurgeon, etc) and

• physiological phenomena (swelling due to osmotic drugs, anesthetics, etc). As a consequence of this intra-operative brain-shift, pre-operative images no longer correspond to reality. Therefore the neuro-navigation system based on those images doesn't necessarily represent the current situation.

In order to face this problem, various teams have proposed to incorporate into existing image-guided neurosurgical systems, a biomechanical modeling to compensate the brain deformations by updating the pre-operative images and planning according to intra-operative brain shape changes. For this, such measured changes (for example the changes of the external shape of the brain tissues in the opening skull region) are given as new boundary conditions to the biomechanical model of the brain tissues that infers the new position of the tumor. Such intra-operative use of a biomechanical model implies that a strong modeling effort must be carried out. Our group tries to provide some clues towards such a modeling effort. Three steps are followed to design the brain model:

• The segmentation of pre-operative images (MRI) to build the external surface mesh of the brain, including the tumor.

• The generation of a volume mesh optimized for real-time simulation.

• The creation of a model of the brain-shift with Finite Elements (FE).

The focus of this paper is the second point. The proposed meshing technique starts from a surface representation of the brain with the tumor and produces a volume mesh of it. Some results about the third point are also shown in section 6.2.

Meshing constraints

Time is crucial in surgery thus, to produce an intra-operative real-time modeling of the brain-shift the FEM computation must be fast. The speed of the FEM directly depends on the number of degrees of freedom the system has, thus an optimal mesh in terms of quantity of nodes must be provided.

A good representation of the tumor as well as the Skull Opening Point (SOP) and the path between them is mandatory because here is where a greater deformation is expected [START_REF] Reinertsen | Validation of vesselbased registration for correction of brain shift[END_REF]. This path will from now on be referred to as the Region of Interest (RoI). A mesh without quality elements can lead to errors in the computation of the FE thus quality is also an important issue in this problem. Therefore the constraints to model the brain-shift in a real-time application are:

Advancing Front

The Advancing Front (AF) technique [START_REF] Lobos | 3d noffset mixed-element mesh generator approach[END_REF][START_REF] Frey | Delaunay tetrahedralization an advancing front approach[END_REF] uses as input a closed surface. All the faces that describe the input domain are treated as fronts and are expanded into the volume in order to achieve a final 3D representation. The selection of points to create the new faces encourages the use of existing points. Figure 1 shows some steps of the AF technique where only a portion of the surface has been taken for clarity purposes. Additional process to improve the quality of the elements can be made. This technique is strongly dependent on the input surface mesh. If the input is accurate, it will have a high number of faces producing a refined volume mesh in the entire input domain. This is contradictory with our first constraint: the mesh must be refined in the RoI and coarse elsewhere. On the other hand, if the input mesh is coarse, this technique should be combined with others strategies of local refinement to respect the first constraint.

In relation to the second constraint, this technique achieves a perfect surface representation because it uses the input surface to produce the final volume mesh. However a little remark can be made: If the input surface mesh is coarse and the RoI is constrained to have a large quantity of points, the resulting final volume mesh will continue to be coarse in the portion of the RoI that intersects the input surface mesh.

Two solutions are proposed: (1) pre-treatment of the input surface mesh to be more refined in the RoI or (2) local refinement in the surface portion of the RoI once the volume mesh is generated.

The third constraint is not managed by this technique. However different algorithms can be applied to improve the quality of the resulting volume mesh, for example in the case of tetrahedral AF meshes, constraint the new tetrahedra to fulfill the Delaunay [START_REF] Delaunay | Sur la sphère vide[END_REF] property as long as the fronts continues their expansion.

Mesh Matching

Mesh matching is an algorithm that starts with a generic volume mesh and tries to match it onto the input domain [START_REF] Couteau | The mesh-matching algorithm: an automatic 3d mesh generator for finite element structures[END_REF][START_REF] Luboz | Orbital and maxillofacial computer aided surgery: Patient-specific finite element models to predict surgical outcomes[END_REF]. The generic volume can be obtained from an interpolation of several sample models. To generate a new mesh, in our case of the brain, the problem is reduced to finding a 3D nonuniform transformation function that will be applied to the entire generic mesh (atlas) and produce in this way the final volume mesh for the current patient. This is shown in figure 2 in a maxillofacial meshing example, from Chabanas et al. [START_REF] Chabanas | Patient specific finite element model of the face soft tissue for computer-assisted maxillofacial surgery[END_REF] Fig. 2. The mesh matching algorithm in a maxillofacial simulation. This technique does not satisfy the first constraint. The atlas does not consider the tumor information thus the mesh is not refined in the RoI. It would be necessary to combine this technique with others strategies of local refinement in order to respect the first constraint. In other words, a generic solution cannot be applied to a problem that is patient dependent, i.e. with positions and sizes of tumor and the SOP that change from one patient to another.

The second constraint is satisfied because the atlas normally resembles the input domain. This can be not true for very specific cases where the brain is malformed.

The third constraint like the second one, is also satisfied because the atlas has a great quality coefficient in the entire mesh. It is a quite perfect model that sometimes has modifications usually made by hand in order to produce the best starting point. A loss of quality can occur only in cases where the input domain and the atlas are not alike.

Regular Octree

The octree meshing technique starts from the bounding box of the surface [START_REF] Schneiders | Octree-based hexahedral mesh generation[END_REF][START_REF]A new approach to octree-based hexahedral meshing[END_REF]. This basic cube or "octant" is split into eight new octants. Each octant is then iteratively split into eight new ones, unless it resides outside the input surface mesh, in which case it is deleted. The algorithm stops when a predefined maximum level of iterations is reached or when a condition of surface approximation is satisfied. The problem with a regular octree mesh for our problem is that it result in a high number of points even in regions where very few deformation of the brain is expected [START_REF] Hu | Intraoperative brain shift prediction using a 3d inhomogeneous patient-specific finite element model[END_REF]. Therefore, a non-optimal mesh would be the input for the FEM producing unnecessary time consumption for the entire simulation. For example the mesh shown in figure 3 is unacceptable for the brain-shift modeling due to excessive quantity of points. Therefore the first constraint is not fulfilled.

The octree by itself does not consider a surface approximation algorithm once the split process is done. Therefore it has to be combined with other techniques in order to achieve a surface representation. Two main approaches are considered:

• Marching cubes [START_REF] Lorensen | Marching cubes: A high resolution 3d surface construction algorithm[END_REF]: this algorithm crops the cubes that lie within the surface and produces, in most cases, tetrahedra.

• Surface projection: this technique projects the points of those elements that intersect the surface, onto it. The main problem is that this can produce degenerate elements unless a minimal displacement is needed.

The quality of the mesh is normally acceptable. The elements with a bad quality can be found only at the surface of the input domain because they have changed in order to achieve the correct representation of the domain.

Delaunay Mesh

In the literature a Delaunay [START_REF] Delaunay | Sur la sphère vide[END_REF] mesh is said to be a mesh of tetrahedra that respects the Delaunay constraint. An example mesh of the Delaunay property can be seen in figure 4 Fig. 4. One output mesh that respect the Delaunay properties. It was generated with TetGen.

A Delaunay mesh can be useful to solve the brain-shift (meshing) problem, however it must be combined with a point insertion strategy that consider the constraints presented in section 2. Current implementations, like TetGen [START_REF] Tetgen | A Quality Tetrahedral Mesh Generator[END_REF], do not offer as much control over point quantity and quality as the constraints needs. One alternative implementation would be the following work flow:

• generate a basic tetrahedralization of the input domain (that doesn't necessarily respect the Delaunay constraint).

• improve the quality of the elements in the RoI and constrain them to a certain size.

• without inserting an excessive number of points, increase the quality of the others elements (outside the RoI).

What is critical in this process is to control the number of nodes. A classical problem of Delaunay meshes is to control the insertion of new nodes when going from one refined zone into another not so refined.

Constraint two can be satisfied by the constrained Delaunay meshes that consider the representation of a domain by using certain predefined points. Those points would be the input surface domain in the brain-shift case. And finally, the third constraint requires the insertion of points to fulfil the underline quality requirement. See also [START_REF] Alliez | Variational tetrahedral meshing[END_REF] for a Delaunay based approach to isotropic tetrahedral meshing.

Discussion Regarding Presented Techniques

The mesh matching and the advancing front techniques are ruled out because it is not easy to incorporate a RoI i.e. to generate a mesh that is refined in a specific region and coarse elsewhere. These techniques need major adjustments in order to achieve the desired mesh.

The Delaunay approach is closer to achieving the desired features, however its emphasis is on the quality of the tetrahedra (i.e. the achievement of the Delaunay property) and not on the quantity of nodes.

Even though the octree technique doesn't accomplish the constraints by itself it allows the refinement of some elements directly by the splitting process. This is a very powerful feature that easily help to fulfill the constraint on local refinement in the RoI. However a process to manage the transition between regions of different refinement levels is needed. A method for surface representation is also needed for the octree in order to achieve an optimal mesh for real-time simulation.

The following section proposes some adaptations provided to the octree technique in order to fulfill all the desired constraints. We consider that this is the best starting point to respect the constraints mentioned in section 2.

Meshing Technique for the Brain Shift

The basic octree algorithm is applied but with two main modifications:

• The condition to stop the refinement is the number of nodes the mesh has.

The quantity of nodes is provided as input to the algorithm.

• One element will be refined only if it is inside the RoI. The first modification corresponds to a basic constraint of real-time simulation. Even with high quality elements, a mesh that would have too many nodes would increase the FE time response during surgery and this is not acceptable.

The second modification is due to the fact that the RoI is the area where the predictions of the developed model should be the most accurate one. Consequently only those elements are refined.

This basic octree mesh is not suitable for the FEM. Figure 5 shows an example of different refinement levels for neighbor elements. The implemented algorithm deals with the different possible cases by adding mixed-elements: pyramids, prisms and tetrahedrons. At the end of this subprocess the mesh is valid for the FEM.

The final step is to achieve a surface representation regarding the quality of the elements. Next part details each step of our algorithm.

Inputs and basic mesh

In order to produce a final mesh several inputs are requested:

• A surface mesh of the input domain (S d ).

• A desired quantity of points. This input will be used as the condition to stop the octree subdivision.

• A surface mesh that describes the RoI (S RoI ). This mesh has to intersect the input domain.

The first step is to obtain the bounding box of the domain. A global search is performed to detect the maximum and minimum coordinates of each axis. The result is an hexahedron that will correspond to the root octant of the octree. Then a point (P out ) that is outside this bounding box is calculated to be used as a reference in several tests that are described later.

The algorithm continues with the classical octree subdivision but with one difference. An octant will be split only if it intersects or is completely inside the RoI (defined by S RoI ). The process continues until the maximum number of nodes, provided as input, is reached.

When an element is created the P out is used to check if the element is completely outside the domain. A virtual segment is created between the P out and each vertex of the element. If one of those segments crosses an odd quantity of S d faces, this element is said to be inside. If not, the element is removed from the mesh. For each new element there is also a test to check if it intersects the S RoI . If not it won't continue to split.

Figure 6 plots the mesh generated for a given quantity of points, the input S d and the S RoI . Note that the mesh is strongly refined in the RoI.

At this step, the mesh is not suitable to carry out a classical Finite Element Analysis (FEA) because of incongruent elements (an element is said to be incongruent when it has one or more points inserted in each face or edge of it). Some modifications of the FE basis functions (see [START_REF] Nesme | Hierarchical multi-resolution finite element model for soft body simulation[END_REF] for example) could enable this mesh to be used for an FEA. However, it was chosen in this paper to remain with a classical FEM approach, which means that the mesh has to be modified in order to manage transitions between incongruent elements. Fig. 6. Octree: only the elements that intersects the RoI are split in order to achieve a higher density in that region: (a) the input and the generated mesh and (b) the generated mesh in solid: dark elements corresponds to the RoI. This mesh has 416 points and cannot be used for classical FEA.

Managing transitions

The developed application has several patterns [START_REF] Hitschfeld | Generation of 3d mixed element meshes using a flexible refinement approach[END_REF] that detect and split the incongruent elements. In particular, we can address to incongruent elements that has at most one point inserted in each face or edge of it. This is a property called "one-irregular". In order to produce a one-irregular mesh, all the octants that do not respect this property are split. In this manner, regions outside the RoI have more points to represent the volume. In other words, the density of points outside the RoI is a consequence of the RoI itself. The resulting mesh can be seen in figure 7.

For each new element generated from this subdivision a test is run to check if it is outside the contour of the brain, as described in the previous section.

It is important to keep the faces as planar as possible and the triangle faces are the only ones that preserve this property all the time. In order to represent the surface the points outside the S d must be projected into it. For this reason, all the elements that intersects the input domain are tetrahedralized. In this manner the outside points can be moved without losing the planar property of the faces.

The octree hierarchical structure is used to assign a default subdivision orientation for each type of element. The neighbor face information is also used to test if there are faces that were already triangulated in which case this orientation is used.

An hexahedron will be split into five tetrahedra as shown in figure 8a. The prism will be divided into three tetrahedra like in figure 8b and the pyramid into only two as figure 8c shows. After the local partition of a surface element into tetraedra, it is possible that some neighbors (such as an hexahedron) get a diagonal in their square face. This diagonal must be incorpored and the neighbor element must be tessellated. This new tessellation depends on the type of the element with a diagonal in their square face. For example, if the element is a pyramid: the base face is triangulated congruently with the neighbors and the propagation (of the triangulation) stops. If the element is a prism it depends on the number of already triangulated faces it has:

• 0: the prism is split in one tetrahedra and one pyramid. The face that wasn't with a diagonal before is updated with the information to be handled by the neighbor element.

• 1: if it is possible to build one pyramid and one tetrahedra this option is preferred and the propagation of diagonalizations stops. If not, one diagonal is added and the corresponding face is updated.

• 2: there is only one option and the propagation stops.

In the case of the hexahedron the inner middle point of it is inserted. With this new point the hexahedron is replaced by 6 pyramids. Then the faces that have diagonals trigger the split of the corresponding pyramid and the propagation stops.

Mesh Quality

An extraordinary number of measures has been proposed, ranging from bounds on solid angles to more complex geometric ratios (See [START_REF] Li | Anisotropic mesh gradation control[END_REF][START_REF] Shewchuk | What is a good linear element? interpolation, conditioning, and quality measures[END_REF][START_REF] Zhang | Adaptive and quality quadrilateral/hexahedral meshing from volumetric imaging data[END_REF] for more details). To our knowledge, only one team [START_REF] Ferrant | Registration of 3d intraoperative mr images of the brain using a finite element biomechanical model[END_REF] dealing with the brain-shift problem through a biomechanical FE modeling of the brain have proposed a mesh quality measure. We decided to use the same measure: the aspect-ratio coefficient (ARC). This coefficient is obtained in different manners for each type of element.

The ARC is normalized so that ARC = 1 corresponds to an ideal element and ARC → ∞ as the element becomes increasingly distorted.

The ARC of the tetrahedra is obtained as follows [START_REF] Ferrant | Registration of 3d intraoperative mr images of the brain using a finite element biomechanical model[END_REF]:

ARC = ( 1 6 6 i=1 (l 2 i )) 3/2 8.47867V el
Where l i (i = 1, . . . , 6) are its edge lengths and V el is the volume of the tetrahedron. The value 1 corresponds in this case to an equilateral tetrahedra.

For the Hexahedron the segments between oppose faces middle points are obtained. The ratio between the longest and the shortest of the three segments will be used as the ARC of this element.

Three segments are also used in the case of the pyramid. The two firsts are constructed using the middle points of opposite edges in the base face. The third segment is the one that represents the height of the pyramid. Like in the case of the hexahedron the coefficient will be the ratio between the longest and shortest segment.

To obtain the ARC of the prism a selection of a square face of it is made. Two segments are obtained as in the case of the pyramid base face. The third segment is obtained using the height of the triangle face that corresponds to the average of the two original triangle faces of the prism. This can be seen in figure 9. The ARC, as in previous cases, is obtained from the ratio between the longest and shortest segment.

Measures over the quality of the mesh are described in section 6.1.

Fig. 9. The aspect-ratio coefficient for the prism is obtained from the ratio between the longest and shortest segment of the three heights presented.

Surface representation

We propose a novel mesh adaptation technique. The goal is to achieve an acceptable surface representation after a quality adaptative octree-based mesh has been generated. The nodes of the mesh obtained until this point of the algorithm do not necessarily lie on the object surface, as described in figure 7.

It is thus necessary to improve the surface representation of the mesh while keeping an acceptable level of quality of its elements. The inputs to our algorithm are:

• The source mesh: a generic volumetric mesh comprising quality elements (in our case the octree-based mixed-element mesh generated above).

• Destination mesh: the surface we want our source mesh to be adapted to (in our case S d .

The key idea of our algorithm is to use a mechanical simulation to constraint the deformation between the source and the destination mesh. We drive the surface nodes of the source mesh with a step-by-step displacement towards the destination mesh, leaving a mechanical model perform the inner nodes relaxation throughout the deformation. It is also important to mention that it was chosen to implement a compressible material behavior (poisson ν = 0.3), in order to allow a mechanical compression of the 3D mesh.

Our approach, unlike direct projection of surface nodes that disregards inner nodes position, is guaranteed to not produce element inversions. To preserve the mesh overall quality we check the elements quality at each step of the deformation. If necessary we artificially increase the mechanical resistance, or stiffness, of those elements that suffer the greatest quality loss before proceeding to the next deformation step. This can result in oscillations between neighboring elements. Therefore in order to guarantee the algorithm termination, a constraint was used: we stop increasing the stiffness of an element if it's Y oung's modulus value reaches a predefined threshold. Note that our goal is not to compute a realistic deformation of the mesh. The virtual mechanical medium is merely used to compute the inner relaxation of the ranges between 1.3 and 2.7, which is quite normal since the octree-based algorithm can generate high quality elements. After the projection, 97% of the elements has a quality value between 1.3 and 7. This can be seen in figure 10. In addition, table 1 shows the quality measurement before and after the projection of the nodes onto the input domain. The average quality values don't drastically change from one to another. The region of interest is mostly represented by hexahedron due to the local refinement of this zone as explained in section 5.2. After projection, the aspect-ratio of the hexahedra is not drastically increased thus the deformation of these elements is reasonable. The same happens with the pyramids and the situation only changes for the tetrahedra, although the aspect-ratio average is still acceptable. We can conclude that the quality in the RoI remains good enough for the simulation. Note that in the presented example no prism was generated. This is because, in our experience, even though they can be found in a mesh, these elements are very rare (just a few patterns consider this type of element).

Modeling

The final mesh is used to simulate the brain-shift with the help of MRI images.

Figure 11 shows how the output mesh is finally put together with the initially scanned images. The tumor is simulated by a sphere and it can be seen how the mesh is congruently more refined in the RoI, i.e. the path between the Opening Skull Point and the tumor.

The FEM deforms the mesh and with this information, the images are updated to show the tissue displacement. A linear elastic small deformation framework is chosen for this computation (E = 1kP a, ν = 0.45). The set of pictures in figure 12 shows how this is done. The first picture shows the initial shape of the brain following the skull aperture. The four points in the picture are control points used to change the model boundary conditions assuming some measurements of the intra-operative changes. The second picture shows the brain deformations induced by such changes in the model boundary conditions, using the control points. And the third picture simulates the resection (removal) of the tumor. As time is crucial in this simulation it is very important to mention that the computation of a deformation is done in 100 ms. This is the process (a) → (b) in figure 12. When there is resection of tissue a main update to the model must be done (re-computation of the elasticity matrix), this operation takes 12 sec for the presented case which is a very acceptable time for surgery. In figure 12 this corresponds to (b) → (c).

Finally figure 13 shows the outputs using different RoIs. In the case of the first one, at the end it has 628 points and the quality coefficient average in the entire mesh is 2.52882. In the second one it has 617 points and a quality coefficient average of 2.9747.

Conclusions

This paper aimed at proposing a method to automatically generate a 3D mesh of the brain adapted to the constraints of an intra-operative use, i.e. a good representation of the area between the targeted tumor and the Opening Skull Point, with a coarser mesh elsewhere in order to allow a fast Finite Element computation. The method was successfully evaluated on a given brain geometry with different simulations for tumor and SOP locations. The algorithm generates a mixed-element mesh that achieves correct surface representation. The overall mesh quality is preserved in the final mesh although the termination rule can lead to excessive quality degradation for some elements. In our experience the degraded elements quantity represents less than 0.1% of the total elements count and a subsequent sliver removal procedure [START_REF] Cheng | Sliver exudation[END_REF] should be applied with little impact on the final mesh topology.

This method could also be used to adapt a mesh in an inter-subject context for fitting a generic mesh on specific patient data. In particular, the issues encountered in the Mesh-Matching procedure as presented in [START_REF] Couteau | The mesh-matching algorithm: an automatic 3d mesh generator for finite element structures[END_REF] and [START_REF] Luboz | Orbital and maxillofacial computer aided surgery: Patient-specific finite element models to predict surgical outcomes[END_REF] should be overcome by the mechanical approach to the displacement of inner nodes.

Fig. 1 .

 1 Fig. 1. The advancing front technique: (a) A portion of a surface mesh with one front to expand, (b) the tetrahedron is created and the new faces can be treated as fronts, (c) another expansion using a recently inserted front and (d) another expansion using already inserted points.

Fig. 3 .

 3 Fig. 3. (a) Input surface triangle mesh of 3152 points and (b) output hexahedral mesh of 20,000 points.

Fig. 5 .

 5 Fig. 5. An example of transition management between different levels of refinement. (a) Two neighbor octants, the left one is more refined. (b) The same two octants but now the right one has a subdivision: 5 pyramids and 4 tetrahedra, making the mesh congruent and suitable for the FEM.

Fig. 7 .

 7 Fig. 7. (a) A one-irregular mesh: each cube has at most one point inserted in each face or edge of it. This mesh has 594 points. (b) A mixed-element mesh that handle the transitions. This mesh has 660 points and can be used for classical FEA.

Fig. 8 .

 8 Fig. 8. (a) The subdivision of an hexahedron into five tetrahedra, (b) The subdivision of a prism into three tetrahedra and (c) The subdivision of a pyramid into two tetrahedra.

Fig. 10 .

 10 Fig. 10. Quality measurement of the elements after the projection.

Fig. 11 .

 11 Fig. 11. MRI images with the final mesh.

Fig. 12 .

 12 Fig. 12. (a) the initial input with skull aperture, (b) brain-shift simulation using 4 control points and (c) update of the model with tumor resection

Fig. 13 .

 13 Fig. 13. (a) a refinement with a RoI in the top-back section of the brain and (b) a refinement with a RoI in the top-middle section of the brain.

Table 1 .

 1 Aspect-ratio per element type average before and after the projection onto the surface.

		Element type Aspect ratio average quantity of elements
		Hexahedra	1.27056	72
	before projection	Prism Pyramid	0 2.1387	0 386
		Tetrahedra 1.58253	1418
		Total	1.68499	1876
		Element type Aspect ratio average quantity of elements
		Hexahedra	1.35187	72
	after projection	Prism Pyramid	0 1.7168	0 386
		Tetrahedra 3.27791	1418
		Total	2.88278	1876

Claudio Lobos, Marek Bucki, Nancy Hitschfeld, and Yohan Payan

Acknowledgement

Francisco Galdames and Fabrice Jaillet for the input surface mesh of the brain used in the explanation of the developed algorithm. This project has been financially supported by FONDECYT 1061227, ALFA IPECA project, FONDEF D04-I-1237 and ECOS-Sud C06E04.

Algorithm 1 Surface representation

Require: source mesh and destination mesh.

Let E be the set of elements in the source mesh, defining the mechanical model. Let S be the set of surface nodes in the source mesh. Let I be the set of inner nodes in the source mesh.

Let D be the destination surface mesh. for all surface nodes P in S do compute the projection vector of P on D: U (P nodes at each deformation step: the Y oungs modulus arbitrary incrementation prevents the most exposed elements from being excessively deformed and possibly degraded.

The entire process can be seen in algorithm 1. Note that in the Do-While loop, if some elements need to be stiffened, the deformation for a given step is redone using the same initial node positions S and I along with the updated elements E. Another very important remark is that this method do not insert new points to the input mesh.

Results

Quality of the mesh

Measurements of the elements quality are performed before and after the projection onto the surface. Before the projection the quality of each element


ANALYSES OF SUBSIDENCE KINETICS ABOVE ROOM AND PILLAR MINES & DEVELOPMENT OF KINETIC CRITERION BASED ON BACK-ANALYSES SUBSIDENCE CASES

The iron ore mines in the Lorraine region of north-eastern France cover a surface area of approximately 3,900 km 2 and involve more than 150 cities. A total extraction method was usually used during exploitation, but this technique was sometimes replaced by the room and abandoned pillar technique when it was necessary to protect surface buildings and infrastructure. Rooms and abandoned pillars were initially designed to ensure long-term mine stability. However, this aim was not always reached and numerous collapses occurred in the Lorraine iron basin, either during production or several years after closure. Several subsidence events were identified during a period of one century; the most recent event occurred in 2012. Two groups of events were identified. The first group represents events characterized by sudden surface movement, and the second group includes event characterized by progressive subsidence.

Numerous studies conducted during the last decade helped improve our understanding the origin mechanism of spontaneous or progressive (slow) events. These studies led to identifying long-term potentially unstable underground zones (hazard zones), particularly deep zones located under buildings and infrastructure. The latest are currently monitored by microseismic networks. A similar high requirement level is needed for unstable zones without taking into account the kinetics of the subsidence, which is very important for safety. Thus, GEODERIS initiated this study, primarily to acquire a better understanding of the kinetics of the progressive events group and identify the main parameters that could explain why some pillars zones developed a maximum amplitude subsidence on surface within a few hours or days, while it took several weeks to months for others. In practice, the results of this study will facilitate adapting the requirement level of microseismic monitoring system to the kinetics of pillar collapse.

To retrieve the real kinetics for each event, historical and sparse documentation was used. Firstly, archives research was conducted, especially testimonies and newspaper articles published during the events. This work helped to gather useful information that sometimes contained a precise description of subsidence duration. The archives research helped to classify some progressive events according to the subsidence kinetics. Secondly, a diagnostic of geotechnical and geomechanical parameters of historic events was realized, combined with a cross analysis of the data. The analysis took into account dimensions of underground workings (e.g., pillar lengths, widths, heights, introduction of a "pillar shape factor", and stress amount). An original expression was highlighted to connect the pillar geometry to the stress state. The kinetics expression was calibrated to classify the historic events taken into account according to the kinetics of movements. The kinetic criterion was applied to predict subsidence kinetics for the potential unstable zones currently monitored by a microseismic system.

KEYWORDS

Iron ore deposit basin, Lorraine region, France, subsidence, kinetic, Roncourt, Auboué, stability, pillars, rooms, stresses.

CONTEXT AND AIM OF RESEARCH

The main iron deposit exploited in Lorraine is from the Aalenian era. This Lorraine oolitic deposit spreads across an area 120 km long and 30 km wide. It surfaces on the hillsides of the Moselle in the east and plunges southwesterly towards the Paris basin, following a slope of 2 to 3°. In the exploited areas, the depth does not exceed 250 m.

In this region, the iron ore was exploited by one or more overlapping layers and several extraction methods were employed. The total extraction method of iron ore was used in areas where there were no buildings, or infrastructure at risk. The almost complete exploitation of the ore was carried out along with the blasting of the mine roof during the works and the surfacing of the subsidence when the area of extraction became wide enough. Hence, subsidence is therefore controlled by this « total » method. Within urban areas the most common method was the abandoned room and pillar extraction method (or partial extraction), therefore preserving housing and infrastructure. This method leaves the pillars in position to maintain the overburden roof. These pillars play an essential role in the long-term stability of the mining structure [START_REF] Chambon | Evaluation de la stabilité d'exploitations par chambres et piliers et comportement du toit. 5 ème Congrès international de mécanique des roches[END_REF][START_REF] Tincelin | La stabilité de la surface à l'aplomb des exploitations souterraines en mines métalliques ou assimilées[END_REF]. In 1997, industrial extraction of the deposit came to an end, leaving almost 40,000 km of tunnels in the Lorraine underground, and a void volume of 1 billion cubic meters.

On numerous occasions during exploitation, the Lorraine Iron Basin (LIB) was at the centre of uncontrolled subsidence having important effects on the surface. The oldest recorded event was the collapse of Audun-le-Tiche in 1902. The most recent one event occurred in 2012 in Neufchef. The history of well-known subsidence phenomena across this basin has brought about the distinction of two subsidence types [START_REF] Didier | Importance of failure mechanisms for management of surface instability risk above abandoned mines[END_REF]El Shayeb, 2001;[START_REF] Tincelin | Effondrement brutaux et généralisés. Coups de toit[END_REF]: those that happen with abrupt surface movement (in a matter of seconds), such as the event in Audun-le-Tiche, and progressive events that happen over a few days, or even months, like the one in Roncourt in 1999 described below. Some of these events (e.g., Auboué in 1996, Moutiers in 1997, and Roncourt in 1999) have been to such a scale that houses have had to have been expropriated because of the extent of the damage.

Over the last decade and in response to the French administration, GEODERIS has carried out many studies to understand the mechanisms for these brutal collapses. A method has been developed and applied to examine the long-term stability of all the former partial-extraction sites in the LIB (Fairhrust, 2003). A retroactive analysis has been carried out on sixteen of the historical events recorded on the LIB to date, half of which are categorized as abrupt, and half as progressive. Geological and geotechnical criteria have been established for the two main types of collapse. For abrupt risk types lying closely to areas with housing and infrastructure on the surface, the risk has been reduced by expropriation. Areas of progressive subsidence are observed by microseismic monitoring systems when the potential damage is substantial. The monitoring program does not consider the kinetics of the phenomena: all the zones are monitored at the same intensity.

SUBSIDENCE KINETICS

Kinetics of Surface Subsidence

Ground subsidence occurs in three main phases (Figure 1). After the initial phase at a slight slow speed, a second phase of acceleration occurs and is followed by a final phase of compaction, during which the speed decreases until it reaches stability. The first effects on the surface, small cracks appearing on the buildings and infrastructure, usually happens just before phase 2; this is the moment where the evolving phenomena can be detected at the surface and hereafter will be referred to as t0. All other time parameters will be in reference to t0. Thus, t1.5 corresponds to the time to reach the critical slope of 1.5%, and tAmax/2 corresponds to the time to reach the amplitude of collapse equal to half of the maximum amplitude expected (or measured). Considering the uncertainty of t0, we defined ν = tAmax/2 -t0 as the time needed to reach Amax/2, and Tc = t1.5 -t0 as the time needed to reach a slope of 1.5%. An event is a danger to public safety if the subsidence causes structural damage to buildings that could affect the security of people. A subsidence is branded as "rapid" if the time separating the first visible signs on the surface and the critical collapse is too short to initiate evacuation. A rapid subsidence must be anticipated by microseismic surveillance of precursors, such as pillar collapse just before surface subsidence. This critical time (Tc) is highly variable depending on which phenomena are being observed on the LIB as illustrated by the following examples.

Roncourt 1999

In the area of Roncourt, damage was reported as soon as November 1998; an increase in reports brought about the establishment of a land-levelling network in February 1999. In 1991 when the commune decided to renovate the road network, a local levelling network was implanted and reference measurement recorded. The latest helped us to assess the subsidence values obtained in February, 1999. So a collapse of In the commune of Auboué, disorder was noticed on the road surface and buildings on 14 October 1996. The levelling network put into place the following day showed a collapse of more than 80 cm in 2 days (i.e., 40 cm/day). Serious damage rapidly affected 57 houses and 70 families, many of whom were evacuated. In this case the collapse was 68 cm, and the Tc was just one day (Table 1, Figure 3). Subsidence measurements for the seven other events included in this study showed that Tc ranged from 1 to 470 days (Table 1). For two events (Doma_01 and Roc_2008), the collapses were sufficiently weak that the damage did not exceed cracks, not impose any danger to people (e.g., Figure 4). Subsidence kinetics is an important element that has an effect on Tc. The subsidence curves can be characterized by a curve-stiffness parameter independent of the final state of collapse. The time ν was chosen as the stiffness characteristic, the advantage being that Tc is less sensitive to the collapse measurements at the beginning of the event. In actuality, after detecting surface movement (damage to buildings), the frequency of measuring levelling markers was increased (generally by 1 measurement per day). This is the date that was chosen as the date of origin from the information gained from witnesses. It is also the origin of the Tc calculation (Figure 1).

The time ν varies from one event to another, similarly to Tc, except for events where the final collapse was too weak (whether it has slow or rapid kinetics). The values of ν enable us to classify the events as either rapid (ν ranges from 1 to 3 days) or slow (ν ranges from 45 to 100 days; Table 1).

Kinetic of Pillars Collapses-Underground Process

The fracture process within mining-works, made up of a group of adjacent pillars, will commence at a point where the security factor is at its weakest, either due to a higher level of stress with the presence of a stacking edge or by a lower level of pillar resistance (brittle pillar fractures). The floor kinetics is the progression of the rupture from one pillar to the next.

The studies conducted in LIB, especially laboratory trials of mechanical resistance to compression, describe the ductile behaviour for this material (Figure 5). Brittle pillar fractures have been dismissed; the rapid kinetics observed for some events is unlikely linked to the mechanism of a brittle fracture [START_REF] Dagallier | Mineralogical and microstructural characterization of the anthropic-origin ageing of iron in Lorraine (France)[END_REF][START_REF] Grgic | A short and long term rheological model to understand the collapses of iron mines in Lorraine (France)[END_REF]. Instead, the collapse process is maintained by the overstress that the adjacent pillars gain; this overstress is proportional to the mesh cutting and the weight of the overlying ground. The level of stress on the panel is therefore an element that is likely to favour the speed of the phenomena. A higher additional solicitation could be behind the quicker propagation speed, enabling the slight local variations of the safety factor to be overcome. The capacity of the pillar to withstand this added stress is linked to its geometry. A more massive and flatter pillar would have a higher compressive strength, expressed by the ratio between the opening of the layer (or the height of pillar) and the width (El1) or length (El2) of the pillar. According to the progression of the collapse, the surface of the pillar can also play a role as much as the latest ratio (El1 or El2) in the opposite direction

Transmission Between the Mine Floor and Surface-Role of Mineworks Overburden

The overburden plays a major role in the transmission of movements between the underground failures and the surface. A theoretical subsidence curve established for the LIB from measurements obtained during exploitation shows surface subsidence in accordance with the collapsed width (D) on the floor (Figure 6). For as long as this width is less than approximately 60% of the depth (H), the effects on the surface are barely noticeable. The shifting will not become serious until the D is more than 80% of H. Hence Tc corresponds to the development of the mine floor fractures between these two lengths. This behaviour is linked to the fairly resistant nature of the overburden of the LIB mines, which contains stiff limestone layers. This phenomenon would be different with a more distortable roof such as the b curve in Figure 6 shows, relating to roofing that has been deformed by the block-caving during a previous extraction of an overlaying stratum [START_REF] Proust | Etude sur les affaissements miniers dans le bassin du Nord et du Pas-de-Calais[END_REF]. The speed of the progression at the front of the pillars collapses will depend on the geometry of the pillars.

RETRO-ANALYSIS OF THE NINE COLLAPSES

The main geotechnical parameters of the nine events are presented in Table 2. Recall from Table 1 that the first five events exhibited "rapid kinetics" and the final four exhibited "slow kinetics". On average, the rapid events were deeper and had wider openings than the slow events. The rapid events are also associated with a regular cutting of identical pillars even if at Crusnes one can see panels with slightly different sizes of pillars leading to different extraction rates (42-53%). However the slow events are more irregular with an irregular mine plan (Roncourt) or a mixture of pillars with different dimensions for the other 3 cases. This result is backed up by the fact that a zone where the pillar dimensions are larger with a lower rate of extraction, yet insufficient to completely stop the progression of the collapse, is capable of stopping for the duration of one month (example of Danté event), 9 months (extension of Angevillers 2009) or even many years as in between Affl_01 and Doma_01, these last two being geographically bordering. The stress within the pillars was usually stronger for rapid collapse (mean 11 MPa) than slow events (mean 9 MPa), which confirms the procedure of weight reporting that is shown above. However, a very rapid event (Mou_1997) had a very low stress level (8 MPa).

Pillar geometry is expressed as the width and length of pillars compared with their height. All pillars had a rectangular shape except for Aub_96_02, where the pillars were diamond shaped. As with the stress, it is not the resistance of the pillar (linked with the ratio El) that plays the leading role in the pillar failure process-given that all the pillars in the area are unstable-but their capacity to take on extra weight during the progression. The ratio El1 calculated on the pillar width does not affect the speed of the phenomenon (see Table 2). However the El2 calculated on the pillar length clearly distinguished the slow events that all had a ratio El2 less than 0,15. Indeed the kinetics of rectangular pillars is slower than that of square ones.

As an attempt to explain the previous observation, we draw in the figure 7 two examples of underground exploitations, one was carried out with square pillars and the other with rectangular pillars; both had the same surface and galleries of equal length. Assuming that a similar procedure had resulted in the fracturing of the first line of pillars and then a second, etc., we can observe that the width of square pillars (D1) at the end of a collapse with the same number of rows of pillars, is greater than the width for long narrower pillars (D2). Assuming that the collapse movements are governed by D, they would require a lower collapsed pillars number with the use of the square. The pillar geometry plays a key-role in the process of fracture kinetics. Therefore we tried to establish an empirical relationship that could account for pillar shape and pillar length. To this end, we introduced two ratio parameters (El1 and El2); the combination of the two was a way to account for the massiveness of the pillars in this case. We recognize that, all other factors being equal, the square-shaped pillars are strongly conductive to rapid kinetics development. An a-dimensional parameter was therefore introduced enabling us to evaluate the variation of any pillar shape with regards to a square shape that its surface would keep. This shape indicator parameter (IF) is between 1 (square pillar) and 2 (an indefinitely long pillar). For any pillar, this ratio is defined as the relationship between the hydraulic diameter (Dh) of the pillar and the length of the pillar (lp).

IF = Dh / lp

(1)

To characterise the massive pillar, we introduced an a-dimensional number made up of the relationship between the pillar surface (Sp) and the exploited opening (w). For this number to be adimensional, we considered the square root of the pillar surface area, which enabled us to obtain the following relationship:

(2)

The formula is no more than the square root of twice the ratio El of the pillar, namely , and .

The implementation of the kinetic indicator (Ic) expressed in equation 2 allowed us to distinguish between the two groups of events. In fact, all rapid events (with a small value of ν) have an Ic less than 6.1 (Figure 8) and slow event have an Ic greater than 6.6 (Table 3).

CONCLUSIONS

Former room and pillar underground exploitations can manifest as subsidence at the surface. When the effects are serious, these events can produce structural damage on buildings. The kinetics of the collapse must then be taken into account. When there is insufficient time between the first signs on the surface (cracks) and structural damage for protective measures to be put into place, surveillance techniques like microseismic recording can be used to predict surface subsidence.

A retro-analysis of nine progressive subsidence collapse events in the last 40 years in the LIB showed that the available reaction time was highly variable, from one day to several months. The stiffness of the collapse evolution curve in relationship to the time was characterized by the time (v) from the beginning of the acceleration of the collapse until half of the maximum subsidence (Amax/2). This length of time is short for the rapid kinetic events and long for the slow kinetic events. Among the nine events, five were categorized as rapid, with v less than 1 week and four events were categorized as slow. The retro-analysis shows that the kinetic is distinguished by the progression of the pillar fracture in the exploited panel. This progression is rapid when the level of stress is high, the panel cutting is regular, and pillars are shaped as close to a square as possible and are short and slim. To quantify the influence of pillar geometry, a shape indicator enabled the separation of the slow and rapid events by a threshold value.

Figure 1 .

 1 Figure 1. Evolution in function of the time taken to reach the amplitude of collapse at a particular point of the Lorraine Iron Basin (LIB)

Figure 2 .

 2 Figure 2. Evolution in function of the time taken to reach the amplitude of collapse measured at the 1026 level in the commune of Roncourt (event Ron_1999) Auboué 1996 "Cité de Coinville"

Figure 3 .

 3 Figure 3. Evolution in function to the time taken to reach the amplitude of collapse measured at marker n 7 in the commune of Auboué (event Aub_96_02)

Figure 4 .

 4 Figure 4. Example of limited damage incurred following the subsidence in Rochonvillers (Roc_2008)

Figure 5 .

 5 Figure 5. Schema of the elasto-plastic ductile mechanical behaviour of Lorraine iron ore

Figure 6 .

 6 Figure 6. Curve schematics of the subsidence evolution according to the collapsed width below LIB: (a) example without extraction underlying the collapsed stratum (b) example with extraction

Figure 7 .

 7 Figure 7. Evolution of the width (D) collapsed at the bottom in the presence of square or rectangular pillars CHOOSING THE KINETIC CRITERIA AND DISCUSSION

Table 1 .

 1 Characteristic of parameters of selected events

	Event	Date	Amax (m)	v (days)	At1,5 (m)	Tc (days)	Kinetics
	Mou_1997	15/05/1997	1.3	1	0.36	1	rapid
	Aub_96_02	14/10/1996	1.2	2	0.51	1	rapid
	Affl_01	30/04/1974	1.35	2	0.77	64	slow
	Doma_01	1976	0.65	3	0.74	N/A	rapid
	Cru_1977	18/11/1977	0.8	2	0.54	26	rapid
	Rue Danté	1972	0.9	45	0.51	41	slow
	Roc_2008	01/01/2009	0.21	170	0.57	N/A	slow
	Ron_1999	01/01/1999	0.75	100	0.42	103	slow
	Ang_2009	4/10/2009	0.68	36	0.53	460	slow

Table 2 .

 2 Geotechnical characteristics of the studied events

	Name of the event	Opening w (m)	Gallery width (m)	Pillar width (m)	Pillar length (m)	ratio (El1/El2)	Extraction rate	H (m)	Total stress (MPa)
	Mou_1997	3.5	6	12	12	0.29/0.29	0.55	120	8
	Aub_96_02	5	6	12	14	0.42/0.36	0.53	170	10.9
	Affl_01	7	6	24	24	0.29/0.29	0.36	255	12
	Doma_01	5	6	27	27	0.19/0.19	0.33	245	11
	Cru_1977	3.8	6	11	25	0.35/0.15	0.50	148	10.9
	Rue Danté	5	8	12	49	0.42/0.10	0.48	170	9.9
	Roc_2008	3	3	10	37	0.30/0.08	0.29	190	8
	Ron_1999	2.5	6	6	85	0.42/0.03	0.53	140	9
	Ang_2009	5	6	11	40	0.45/0.13	0.44	175	9.4

Table 3 .

 3 Characteristics of the nine events chosen for kinetic analysis Projection of the 9 events on the plan (Ic, stress) (red: rapid; blue: slow)

	Event	IF	Ic	Kinetics (with time ν)
	Mou_1997	1	3.43	Rapid
	Aub_96_02 1.08 2.79	Rapid
	Affl_01	1	3.43	Rapid
	Doma_01	1	5.4	Rapid
	Cru_1977 1.39 6.06

Rapid Rue Danté 1.58 7.79 Slow Roc_2008 1.57 10.1 Slow Ron_1999 1.87 16.88 Slow Ang_2009 1.57 6.58 Slow Figure 8.

cm had already occurred at a speed of approximately 0.3 cm/day. Movement continued at a maximum speed of 2 cm/day. This very slow subsidence caused structural damage to houses that required the evacuation of buildings during 1999. In this case the critical collapse was 42 cm (with a slope at 1.5%) and the critical time (Tc) was 103 days, the event can therefore be characterized as slow (Table1, Figure2).


Long distance propagation of shock waves in the open atmosphere

C. Proust, L. Heudier, INERIS, Verneuil-en-Halatte, France K. Vilalta, Nobelclad Europe, Rivesaltes, France

Context and Objectives

Despite significant progresses in the understanding of the formation and propagation of pressure waves in the open atmosphere in the past decades, experience form past accidents (AZF, Toulouse, 2001, Billy-Berclau, 2003 and even Buncefied 2005) reveals that overpressure levels may be greatly ill estimated. This is obviously a concern when the protection of citizens is concerned.

There are some indications that both the topography and weather conditions may play a significant role but this remains rather qualitative and there is a need to clarify such points. In this paper, the results of a preliminary experimental campaign are presented in which up to 500 kg eq TNT was detonated in the open atmosphere and the pressure field measured up to 10 kms.

Existing knowledge

Most traditional models for the propagation of shock waves in the open atmosphere assume an isotropic and homogeneous atmosphere ([1], [START_REF] Van Den Berg | The Multi-Energy Method -a framework for vapour cloud explosion blast prediction[END_REF]) but their results may diverge significantly especially when the overpressure is small (Figure 1). When analyzing the damages following an industrial ( [START_REF] Branka | Intervention de L'INERIS après le sinistre survenu le 27 mars 2003 sur le site de BILLY-BERCLAU de la Société NITROCHIMIE[END_REF][START_REF] Mouilleau | Initial analysis of the damage observed in Toulouse after the accident that occured on 21st of September on the AZF site of the Grande Paroisse company[END_REF]), it appears that the evolution of the pressure field at large distances seem to be different from the traditional 1/x decay with the distance. In fact (Figure 2), the decrease of the overpressure seems much less especially in the range 0-50 mbar. In practice from one method to the other a ratio of 2 at least are possible especially in the low pressure range (below 100 mbar).

Several reasons may explain these discrepancies. The vertical gradient of the wind may serve as a "wave guide" which may explain the different-from-isotropic propagation of the shock wave [START_REF] Baker | Explosions in air[END_REF]. Second, density gradients may also play a role by diffracting the waves by this point does not seem to have been deserved a great attention (Figure 3 from [START_REF] Rybnov | Experimental studies of the effect of the ground atmospheric layer and the underlying surface on the amplitude of weak blast waves from ground chemical explosions[END_REF]). The topography and the nature of the ground is potentially also of some importance for instance by modifying the frequency spectrum of the wave (low frequency waves are less absorbed by the atmosphere).

Experimental site and instrumentation

The experimental site is belonging to Nobelclad Europe and is located in the south of France, near Perpignan. This site is a "pyrotechnical" zone devoted to the everyday activity of the company specialized in the welding of thick metal plates using the heat and pressure of a detonating material (Figure 4). The explosive is a powdered ANFO layered in typical thicknesses of 90 mm. Up to 1000 kg can be detonated corresponding to an energy of about 500 kg TNT [START_REF] Nobili | The Explosion Bonding Process[END_REF].

Results

A typical example is presented below corresponding to the detonation of 333 kg of TNT layered over 3 plates (thickness of the ANFO layer=95 mm, test n°00052).

By comparing both the shape of the pressure traces and the characteristic parameters to the theoretical estimations (in figure 7 from []), it is clear that the actual shocks waves are significantly different. The overpressure about 30 m (K1) from the detonation initiation is lower that the theory but the duration and impulse are in line. This suggests that most of the initial energy of the ANFO charge is transferred to the wave. The shape of the wave is not "triangular" suggesting the wave is not totally formed or was modified by the local topography. At larger distance (K2@300 m), the overpressure and impulse are greater than predicted but this may result from the channeling effect of the "corridor". Very far from the explosion platform (figure 8), the shape of the waves is not any more a shock. The initial shock wave has lost the high frequency portions (max frequency at 1400 m is about 100 Hz), presumably absorbed by the air humidity [START_REF] Kinsler | Fundamentals of Acoustics[END_REF]. Note that the rarefaction wave following the initial shock and the subsequent "vibrations"/echoes do not attenuate as rapidly. There is apparently a significant directional effect (compare 3800 and 4200 m) although there was nearly no wind during the test. Lastly, the maximum overpressure seems to be much below the acoustic approximation [START_REF] Koper | Empirical scaling laws for truck bomb explosions based on seismic and acoustic data[END_REF].

@1400 m @3800 m @4200 m @9400 m Some evolutions of the overpressures as function of the explosive charge are presented on figures 9 (near field) and 10 (far field). A considerable scattering of the results is seen which is greater, the larger the distance. From the available data, correlating the far field overpressure measurements to the global weather parameters does not provide any clear relationship.

Conclusions

In this work, an effort is made to clarify the environmental parameters on the propagation of pressure waves at larges distances (kms). Experiments are currently being performed on a pyrotechnical site detonating hundreds of kgs of ANFO. Preliminary results are presented but more data will be issued shortly.

Present information reveals a very large scattering of the data and that the theoretical models are not in line with the measurements.

Figure 1 :

 1 Figure 1 : different models for the spherical propagation of a shock wave produced by detonating 1000 kg TNT. (numeric : Multi-energy method ; experimental : TM5-1300, French regulation).évolution de la surpression estimée en fonction de la distance à l'épicentre

Figure 2 :

 2 Figure 2 : upper (squares) and lower (diamonds) bounds of the pressure levels (from the damages) as function of the distance after the Billy-Berclau explosion (50 kg of explosive)

Figure 4 :

 4 Figure 4: explosive cladding of two metal plates

Figure 6 :

 6 Figure 6 : Far field instrumentation

Figure 7 :

 7 Figure 7 : overpressures measured in the near field (330 kg TNT, test n° 00052)

Figure 8 :

 8 Figure 8 : overpressures measured in the far field (330 kg TNT; test n° 00052)

Figure 9 :

 9 Figure 9 : nearfield overpressures as function of the reduced distance


Ben Moussa

Fortin

Thermal Radiation Contribution to Metal Dust Explosions

Introduction

Dust explosions remain a safety and scientific challenge. There has been a gradual evolution in the prevention and mitigation of dust explosions over the past twenty years [START_REF] Eckhoff | Current status and expected future trends in dust explosion research[END_REF][START_REF] Amyotte | Application of inherent safety principles to dust explosion prevention and mitigation[END_REF]. This evolution results partly from a better understanding of the phenomena involved; namely the properties and characteristics of dust explosions [START_REF] Cashdollar | Overview of dust explosibility characteristics[END_REF], the propagation of the flame through dust clouds [START_REF] Proust | A few fundamental aspects about ignition and flame propagation in dust clouds[END_REF], the ignition and the combustion time of particles [START_REF] Huang | Combustion of bimodal nano/microsized aluminum particle in air[END_REF][START_REF] Goroshin | Emission spectroscopy of flame fronts in aluminum suspensions[END_REF][START_REF] Yetter | Metal particle combustion and nanotechnology[END_REF] and the incidence of the size and concentration of particles [START_REF] Di Benedetto | Modelling the effect of particle size on dust explosions[END_REF]. Today, the modeling of this type of explosions is generally derived from the modeling of gas explosion because of the similarities between the flame propagation processes in both media at least for some categories of dusts. Flour, starches and sulfur powder certainly belong to this panel because it was shown that the particles are gasified ahead of the combustion zone and that heat conduction through the front leads to this transformation [START_REF] Proust | A few fundamental aspects about ignition and flame propagation in dust clouds[END_REF]. The situation might be different if the particles do not evaporate and/or if the major part of the combustion process is heterogeneous. This may occur with fine metal particles as aluminum. In particular, the combustion process might be heterogeneous and results in the presence of solid residues at very high temperature [START_REF] Marion | Studies on the Ignition and Burning of Levitates Aluminum Particles[END_REF] potentially transferring most of the thermal energy to the reactants (aluminum particles) by thermal radiation [START_REF] Blouquin | Contribution à l'étude théorique des interactions entre combustion et rayonnement[END_REF].

Only very few experimental observations and theoretical considerations are available but seem to confirm that thermal radiation exchange is significant in the thermal balance of the flame [START_REF] Cassel | Factors affecting flame propagation through dust clouds[END_REF] and that it could lead to a dramatic acceleration of the flame [START_REF] Proust | Amplification des effets des explosions sous l'effet du rayonnement thermique[END_REF]. The aim of this work is to investigate further this phenomenology.

Theoretical Approach (numerics)

The coupling between thermal radiation, thermal conduction and particle combustion is extremely complicated and challenges the analytical approach if no reasonable simplifications can be made. In the absence of sufficient convincing experimental results, it is difficult to sort out the potential assumptions. And without some reasonable theoretical grounds, it is difficult to orientate the experimental work. It was decided to investigate the theoretical aspects in the most flexible and open manner. The numerical resources of the Discrete Element Method [START_REF] Fortin | Numerical simulation of granular materials by an improved discrete element method[END_REF] were used for this purpose under the frame of MULTICOR software developed by LTI [START_REF] Fortin | Simulation numérique de la dynamique des systèmes multi-corps appliquée aux milieux granulaires[END_REF].

It was attempted to simulate the radiative transfer between particles in suspension by reproducing a simple total absorption model (without reemission) of the flux impacting each particle. It's schematically the Beer-Lambert attenuation low used by Essenhigh [START_REF] Essenhigh | The thermal radiation theory for plane flame propagation in coal dust clouds[END_REF] and by many authors after him. The details of this one-dimensional attenuation law (under the assumptions considered) are given in the works of Essenhigh [START_REF] Essenhigh | The thermal radiation theory for plane flame propagation in coal dust clouds[END_REF].

Every ignited particle is considered a point source that radiates with a power, Q Ri , equivalent to the power of the entire flame:

4 f i s i Ri T S Q   
where ɛ i is the emissivity of the particle Ω i , σ s is the Stefan-Boltzmann constant, S i is the surface of the particle Ω i and T f is the flame temperature (≈ 3000 °k).

Based on this model, the radiative heat flux density, Q ij , received by the target particle varies inversely to the square of the distance separating the source and the target as follows:

2 . 4 ij Ri ij d Q Q  
where d ij is the distance separating the particle source Ω i and the receiving particle Ω j as illustrated on Fig. 1. We take into consideration, in addition to the remoteness of particles from each other, the attenuation of the radiation absorbed by the screening effect of the particles between transmitter and receiver. The heat flux density received by the target may be then expressed as follows:

                                    1 1 2 1 2 2 2 2 1 1 ) ( 4 1 . ) .( 4
S d Q d Q S Q Q    
where D is the mean distance between two particles.

The proposed results are obtained after programming this model with the DEM code MULTICOR [START_REF] Fortin | Simulation numérique de la dynamique des systèmes multi-corps appliquée aux milieux granulaires[END_REF]. Dust particles are considered as discrete elements. Fig. 2 details the 2D simulation of a uniform dust cloud of 121 spherical particles (d p = 10 μm) distant of 100 μm from each other. The dust concentration is equal to 600 g/m 3 . Air conductivity is considered equal to zero so that the only heat exchange mode taken into account is radiation. Particles located at the plane (x = y = 0) are initially considered as heat sources at the flame temperature (≈ 3000 °K). We suppose that as soon as the aluminum melting temperature (≈ 933 °K) is reached, particles ignite and their temperatures rise instantly to the flame temperature. As shown in

Experimental Work

In addition to the numerical approach, an experimental work is being conducted. From a simple analysis [START_REF] Proust | Amplification des effets des explosions sous l'effet du rayonnement thermique[END_REF], and waiting for more quantitative results, the interesting range of particle diameters is 1-10 μm. The following particles were selected (table 1).

In these size ranges, the interparticle cohesive forces are very strong and agglomerates are formed. The various techniques of disaggregation were studied knowing that the suspension should be laminar in

Q ij [W]
x

[mm]

Radiative flux measurement y x respect of the simulation conditions. The study confirmed that the fluidized technique of a bed of particles with elutriation is the most promising.

Conclusion and perspectives

The modeling of radiative heat exchange trough the flame front is particularly complex because in order to be complete, it must integrate light scattering by particles (reflection, transmission, diffraction), absorption and emission. The coefficients of absorption and emission are dependent on both the wavelength and particle size. A local analysis of the radiation transfer process is still needed to calculate these coefficients but also to determine the incorporation of the diffusion through the flame front. Experimentally, succeeding in dispersing fine aluminum particles will enable the measuring of the burning rate, the particle temperature and the heat radiation in respect of the simulation conditions. We will try to create the theoretically predicted conditions for instabilities induced by radiation.

Figure 1 .

 1 Figure 1. Radiative exchange between particles in suspension

Fig 2 ,Figure 2 .Figure 3 .

 223 Figure 2. (a) Dust temperature and (b) temperature contours at t = 0.5 s

  is shown on Fig.4. It basically consists in a fluidized bed device[START_REF] Goroshin | Emission spectroscopy of flame fronts in aluminum suspensions[END_REF] feeding the flame propagation tube (3) by the bottom while the top end is open and communicates with a separator (1) via a valve[START_REF] Amyotte | Application of inherent safety principles to dust explosion prevention and mitigation[END_REF]. The ignition device is located at 15 cm from the base of the tube (4). The dust-air mixture is fed into the tube during a few tens of seconds. After this sequence, within a few seconds, the air flow is stopped, the top valve is closed and the fluidized bed device is removed and replaced by a flexible tube communicating with the exterior and the igniter is triggered. The flame propagates freely from the bottom to the top giving access to the burning velocity. A typical example is shown on Fig.4. A number of "thermal" measuring techniques were tested and the optical measurement of the monochromatic heat flux seems to be the best solution (again an example of spectrum is shown on Fig.4).

Figure 4 .

 4 Figure 4. Flame propagation tube and instrumentation

Table 1 :

 1 Provided data about the nature of aluminum powders

Acknowledgements

This study is carried out under RADIAN project (2013-2016). The authors would like to thank the "Région de Picardie" for its financial support.


M Al Heib

C Duval

J.-M Watelet

P Gombert

Back-analysis of the collapse of the Clamart chalk underground quarry -Paris (France) Rétro-analyse de l'effondrement de la carrière de Clamart -Paris (France)

Many large collapses have occurred during or after underground quarry operations. The collapse of the Clamart underground quarry (in 1961, South of Paris-France) is well known as one of the most catastrophic. This paper presents a back-analysis of the case study using an empirical approach. Based on a set of assumptions, the empirical approach has indicated a reliable scenario that may explain the collapse: the large collapse may have been caused by a lack of bearing capacity of the pillars and the sudden failure of a hard limestone bed in the roof, the load on which would have been concentrated until overload occurred. The bed weakness between the two quarrying levels may also have had an influence on the general lack of stability of the quarry. Moreover, under the given hypothesis, the study shows a very low influence of the water table due to flooding of the underground quarry. RÉSUMÉ Les effondrements peuvent affecter les cavités souterraines pendant et après le creusement. L'effondrement historique en 1961 de la carrière de Clamart (sud de Paris -France) est connu comme le plus catastrophique. Cet article présente une rétro-analyse basée sur une approche analytique-empirique. Les résultats ont montré que la cause principale de l'effondrement est l'insuffisance de la résistance des piliers et la rupture brutale du banc de calcaire du haut toit. Les piliers ont été stables grâce à la présence du banc de calcaire dans le toit. D'autres facteurs comme l'ennoyage et la rupture de l'intercalaire ont également contribué à l'effondrement, avec un rôle moins important.

INTRODUCTION AND OBJECTIVES

Man-made and natural shallow cavities have caused many different problems with geohazard, and sometimes large collapses. In France and around the world, the old quarries were often excavated using the room and pillar method (Figure 1). This method results in substantial geohazard long after the underground quarry has been abandoned and therefore poses a serious public safety issue today [START_REF] Gutiérrez | Investigation of a large collapse sinkhole affecting a multi-storey building by means of geophysics and the trenching technique (Zaragoza city, NE Spain)[END_REF]. Shallow man-made chalk cavities (depth < 50 m) can be affected both by local collapses and large collapses (Vachat, 1982[START_REF] Kaufmann | Geohazard map of cover-collapse sinkholes in the 'Tournaisis' area, southern Belgium[END_REF], Gilbride et al., 2005[START_REF] Edmonds | Karst and mining geohazard with particular references to the chalk outcrop, England[END_REF]. The impact of collapse on the surface can be subsidence or sinkholes. Such collapses can be found in the literature [START_REF] Bell | Unusual cases of mining subsidence from Great Britain, Germany and Colombia[END_REF]. A large collapse covers a large surface and can culminate in a small earthquake or rockburst; in one gypsum underground quarry, the local magnitude was registered at near 3.1 [START_REF] Wang | Investigation of catastrophic ground collapse in Xingtai gypsum mines in China[END_REF].

In France, there are about 150,000 known quarries, and many of them are located in the chalk, limestone and gypsum of the sedimentary Paris Basin [START_REF] Ph | In-situ and laboratory tests to evaluate the impact of water table fluctuations on the stability of chalk underground quarries[END_REF].

Chalk is a white and soft limestone which forms the base of the Parisian sedimentary basin. The chalk was extensively underground quarry until the beginning of the 20 th century. The room and pillar mining method was used to protect the surface of the land and to ensure the long-term stability of the pillars. However, the long-term stability is never guarantee, and 391 subsidence events related to underground man-made cavities were recorded from 1977 to 2000. face more than 500 m 2 ) from 1827 to 1991.

SUDDEN GROUND COLLAPSE AND SUB-SIDENCE

The disastrous collapse of the chalk underground quarry of Clamart and Issy-les-Moulineaux, two bordering suburbs of Paris, on June 1, 1961 is one catastrophic example (Figure 2). The collapse of the Clamart underground quarry, where mining ceased in the 1880s, caused the abrupt subsidence of three hectares of ground surface and the destruction of a whole urban district, including residential buildings, several roads and sports facilities (Figure 3). Twenty-one people died and more than fifty people were injured.

According to seismographic reports, the ground collapsed in two stages. The southern part, underneath the football stadium (Figure 3), collapsed first, followed by the northern part about half an hour later. The centre of the southern zone suffered the greatest subsidence, at 6 metres. The centre of the northern zone has dropped by 4 metres. Large cracks of about 1 metre wide developed on the edges of the area of the collapse and between the southern and northern zones. As a result of the subsidence, some buildings in the area suffered major damage or were totally destroyed. All buildings shown in black on the figure had to be demolished for safety reasons.

OVERVIEW OF MINING

After having been operated for more than 100 years, the Clamart underground quarry was abandoned before 1880 due to the reserve limitation and mining licence, and the exact date is unknown. In the collapse zone, two extraction chalk levels were operated, both of which employed a room and pillar method. Different extraction ratios were followed; the galleries on the upper level were about 7 meters high in both the northern and southern zones. The pillars between the two zones were irregular in shape, showing a hybrid design (Figure 4). The rooms on the lower level were about 4 metres high and existed only in the southern part, as shown in Figure 5. On the two levels, the area between the northern and southern zones was probably situated on a property line. The interlayer bed (immediate roof of the lower level) was about 4 metres high, and the two levels were perfectly superimposed. The geometrical characteristics of the different underground quarry areas shown in Figure 5 can be found in Table 1. The extraction rate ranges from 36% to 68%. The pillar width/height ratio varies between 0.75 and 1.4. Irreg.

-4*5 36

3.1 Objective Various experts have attempted to explain the collapse between 1961 and 1979 on the basis of in-situ observations. The causes leading to the collapse have not yet been fully explained, although several hypotheses have been proposed in reports and articles dealing with this disaster. IGC considered the main cause of the collapse was the failure of the limestone bed located on the upper roof (unpublished expert reports from 1961, 1976). Maury [START_REF] Maury | Effondrements spontanés. Synthèse d'observations et possibilité de mécanisme initiateur par mise en charge hydraulique [Summary of observations and possibility of initiator mechanism by hydraulic loading[END_REF] considered the main factor of the collapse to be the hydraulic load and water pressure induced by rain many days before the collapse. The aim of this article is to allow a better understanding of the origin of the collapse through an evaluation of the stability of the underground quarry and an analysis of triggering factors and their associated failure mechanisms.

ASSESSMENT OF UNDERGROUND QUARRY INSTABILITY

From the date of the collapse up to now, several communications by experts and authorities attempted to explain the collapse of the Clamart underground quarry [START_REF] Maury | Effondrements spontanés. Synthèse d'observations et possibilité de mécanisme initiateur par mise en charge hydraulique [Summary of observations and possibility of initiator mechanism by hydraulic loading[END_REF]. They considered the explosion, slipping of the limestone layer on the plastic clays, weak pillars under the water level, hydraulic pressure in the cover layers after previous pipe leaks or heavy rain. To this day, it remains difficult to determine the real cause of the instability. In order to achieve the objective, the adopted methodology is divided into two parts; the first part concerns the stability of the underground quarry before the collapse, while the second part concerns the triggering factors.

An empirical-analytical approach was used in parallel to identify the relevant causes of the catastrophic collapse.

In order to assess the stability, a main hypothesis is made: the large collapse of the underground quarry is supposed to have been initiated by the failure of the pillars in the area where there is the deepest subsidence and where the extraction ratio is near the maximum. This corresponds to the area in which the two levels are superposed. The collapse of the pillars could be induced by secondary factors: collapse of the immediate roof, hard limestone bed and water, and underground quarry flooding.

Assessment of pillar stress

The analytical-empirical approach allows for the successful design and the stability assessment of the underground quarry, depending on the individual pillar stability [START_REF] Bekendam | Stability and Subsidence Assessment over Shallow Abandoned Room and Pillar Limestone Underground quarries, Engineering Geology for Infrastructure Planning in Europe[END_REF]. The average pillar stress can be estimated by the tributary area method based on the following parameters for the geometry of rooms and pillars, the depth and the specific weight of overlying rocks. The method was applied to the Clamart zones (Table 2). The minimum pillar stress calculated out of the 6 zones of interest is 1.75 MPa (Zone ZSS1 on Figure 5). This corresponds to the area of maximum subsidence, thus reinforcing the initial hypothesis of the initiation of the collapse.

Pillar failure occurs when the pillar is loaded beyond its bearing capacity and the safety factor is less than 1. A large-scale collapse can occur as a catastrophic collapse when an excessive convergence exists; the slenderer the pillar is, the wider the area of collapse is (Esterhuizen et al., 2006). The pillar bearing capacity was first assumed to be the same as the chalk's unconfined compression strength (UCS). This was possible thanks to the compression testing that had been conducted (with different values of water content) just after the collapse (in June 1961) and resulting in Table 3. We adopted the weakest value of compression strength (2.1 MPa) for the collapse analysis.

A simple comparison between the unconfined compression strength and the pillar stress indicates that the pillars should not have collapsed, except for ZNS1 (factor of safety <1). However, the subsidence of the surface above the ZNS1 was less than in other zones. In addition, it is known that laboratory strength does not take the large-scale effect into account, so the pillar strength is in fact lower than the laboratory strength. The compression strength of the chalk on a large scale has been estimated at 1.9 MPa based on laboratory tests and back-analysis on the region's chalk underground quarries in the Paris region [START_REF] Watelet | Méthode d'analyse et diagnostic des conditions de stabilités de carrières souterraines [Method for analysing and diagnosing stability conditions of subterranean quarries[END_REF]. And yet it is still too high to explain the collapse. According to our observations and the known relationship between the underground quarry opening and the subsidence, the maximum subsidence was located in ZSI1 (Figure 5); one can assume that the collapse of the pillars started there. In order to have an idea of the pillars' strength, we assume the factor of safety (FoS) of zone ZSI1 (lower level) to be 1. It must be noted that the safety factor may have been lower (meaning a lower strength of the pillars), but the underground quarry stabilized temporarily for an unknown reason (the resistant limestone roof and/or the local confining pillars by backfilling). Anyhow, given that the pillars were doubtlessly unstable (since they eventually collapsed), the safety factor is 1 at the most and the ultimate strength of the pillars is 1.75 MPa. Considering this value as the failure criterion:  the failure criterion is exceeded in zones ZSS1, ZSS2, ZSI2, all of which also collapsed;  the failure criterion is not reached in zone ZSI3, which is still standing;  the failure criterion is far exceeded in zone ZNS. This zone has not collapsed, probably because of the reinforcement on the pillars observed in this area. This probably means that they were showing signs of instability and consequently reinforced.

However, it must be noted that, before the collapse, the upper level showed no disorders due to overloading (no fissures and cracks) and had never been treated by any reinforcement (IGC, 1961). This is inconsistent with the instability of the pillars of ZNS1 and ZNS2, unless no existing disorders were observed in the chalk. Indeed, signs of instability are known to be difficult to see in chalk. It follows logically, however, that they should have been noticed the same way that damage was noticed in area ZNS.

Hydrogeological context of the 1961 event

May water be considered as one instability factor that led to the collapse? The underground quarry collapse occurred when the chalk water table was at an average or moderately high level, following 15 years of low water table. The lower underground quarry level had probably already been flooded for a long time after the abandonment of the Clamart underground quarry (1880).

Concerning the hydrological context, the level of the Seine can vary from a minimum of around 26 m a.s.l. (inducing chalk drainage) to a maximum of 35 m a.s.l. (inducing chalk recharge). The Seine level might have been low (around 27 m a.s.l.) during the months preceding the event.

In terms of rainfall, 1961 was a dry year (540 mm) following a rainy one (>700 mm in 1958 and 1960). Looking at a monthly scale, May 1961 is one of the rainiest months of the 1941-1960 period. It is mainly characterized by heavy rainfalls concentrated on a small number of days. Other potentially aggravating factors related to water were recorded shortly before the collapse:  earthworks were underway in eastern edge of the football stadium and caused the rupture of a water pipe, releasing about 100 m 3 of water above the mining area within 2 hours;  mishandling occurred at the drain valve sprinkler system of the stadium, located in the northwest end; this resulted in the spillage of 2500 m 3 of water during the first 3 weeks of March. The global hydrogeology context of the 1961 event can be summarized as follows:  the chalk aquifer water table was average since the beginning of the year;  the Seine had already begun its summer low flow and low level period;  the rainfall was characterized by a dry first semester followed by a very rainy last week of May, the week preceding the event.

CONCLUSION

The main objective of this back-analysis of the large collapse of the Clamart chalk underground quarry was to identify and to understand the origins of the tragic event in order to predict such events for identical configurations. The back-analysis indicates that the chalk underground quarry was, before the collapse, very close to its ultimate strength and collapse. The empirical-analytical approach indicates that the strength of chalk must be more than 1.75 MPa (the pillar stress) for stable pillars; the underground quarry may be unstable because the laboratory tests and site effect indicate that the strength is very close to the actual pillar stress.

Therefore, the scenario of the large underground quarry collapse probably occurred in the upper pillars. However, no damage was observed on the upper pillars despite being close to the ultimate load. This puts forward the role of the resistant limestone bed (Figure 4): it would have taken back the overload, preventing pillar damage, and it then would have lost its stiffness, leading to the overloading of the pillars and the collapse of the underground quarry. The causes of the limestone bed failure may be multiple origins: hydraulic overload fracture, progressive deterioration of the bonds to the minimal recovery due to the water, and local failure of a pillar sufficiently increasing the overload on the roof.

The bed between the two chalk levels may have fallen during the general collapse, explaining an observed collapse on two levels and affecting the strength of the pillars.

The scenario explaining the breakdown of the upper pillars only concerns the fall and the resumption of stresses through the roof. It does not examine the failure mechanisms in the roof. The presence of two exploitation levels has some importance, but the main role concerns the weakness of the chalk layer and the existence of the hard bed in the overburden.

The methodology developed in this paper uses a scientific and operational approach to explain the causes of such phenomena and to identify the role of different factors suspected by experts. The backanalysis and the knowledge of local context put together are very useful for the understanding of the complex geotechnical phenomena.

Figure 1 .

 1 Figure 1. Room and pillar mining method for shallow underground quarry. The majority of collapses concerns gypsum underground quarries with a surface subsidence of less than 2000 m 2 . The chalk underground quarries around Paris give cause for concern of large-scale collapses: Château-Landon (7 deaths in 1910), the Clamart collapse (Figure 2) -largest one at 30 000 m 2 (21 deaths in 1961), Chanteloup (1 death in 1991) and Bagnolet (2 deaths in 1993).

Figure 2 .

 2 Figure 2. Location of underground quarries and karsts (in Paris).

Figure 3 .

 3 Figure 3. The districts of Clamart and Issy-les-Moulineaux after the ground collapse (IGN, August 1961).

Figure 4 .

 4 Figure 4. Simplified cross-section of the overburden chalk underground quarry.

Figure 5 .

 5 Figure 5. The upper and lower levels of the underground quarry (left and right, respectively, from IGC report).

  N N

Table 1 .

 1 Geometrical characteristics of the upper and lower levels.

	Zone	Pillar	Pillar di-	Gallery	Extraction
		shape	mensions	dimension	ratio %
			(l*L m)	(W*H m)
	ZSS1	Rect.	7*9-15	7	57
	ZSS2	Hybrid	5*5	4*7	55
	ZNS1	Square	7*23-27	4*6.8	68
	ZSI1	Rect.	7*18	4*5	45
	ZSI2	Rect.	7*14	4*5	51
	ZSI3

Table 2 .

 2 Pillar vertical stress calculated by the tributary area method for each zone from Figure5(safety factor was calculated for unconfined compression strength equal to 2.1 MPa)

	Zone	 (%)	D (m)	 (t/m 3 )	p(MPa)	F
	ZSS1	57	42		1.95	1.07
	ZSS2	55	42		1.87	1.12
	ZSI1	45	52	2	1.75	1.20
	ZSI2	51	52		1.96	1.07
	ZSI3	36	52		1.50	1.40
	ZNSI	68	35		2.19	0.96
	: extraction ratio, D: depth, : density, p: pillar stress, F: safety
	factor

Table 3 .

 3 Unconfined compression strength (UCS) characteristics of the Clamart chalk underground quarry(IGC, 1961, unpublished) for different values of water content (W)

	W	UCS (MPa)	UCS av. (MPa)
	50%	5.6; 3.5; 3.4	4.17
	65%	5.0; 2.8; 2.5	3.43
	83%	3.2; 2.0; 2.4;	2.53
	90%	2.7; 2.5; 2.0	2.4
	Immersion for 3	2.1	2.1
	weeks