Francis Wojtkowiak
email: francis.wojtkowiak@ineris.fr

Jean-Christophe Couillet

Xavier Daupley

Geotechnical and environmental impacts on the surface of the water rising in French Underground coal mines after closure

The water rising in Underground coal mines, after their closure and, therefore, to the interruption of the dewatering pumping, may have non-negligible impacts of various types on the surface and on the environment. In the scope of numerous studies conducted in the various French coal basins, INERIS developed a study method to define the residual mining risks, äs well äs their potential impacts on the environment, resulting from the definitive closure of the coal mines. Amongst the main consequences of the flooding of mines on the environment, we will more especially deal with the impacts of the water rising on the stability conditions of the surface (risk of subsidence restoration or elevation of the ground; risk of collapse ...) and with the risks of gas emissions to the surface or even accumulations in places not or badly ventilated. The various phenomena and mechanisms that are at stake and the compensatory measures the most currently implemented to reduce and/or manage these risks are briefly depicted. Finally, a mathematical model to efficiently manage mine gas resource after the mine closure, while guaranteeing the safety of people and facilities and the protection of the environment, is more thoroughiy detailed.

INTRODUCTION

The definitive closure of an Underground mine generally gives rise to the interruption of the dewatering pumping. The resulting rise of the water level within the mine may have non-negligible effects of various types on the surface. The major effects and associated risks are äs follows: 1. Potential impacts on the quality of Underground and surface waters and the appearance of water outlets more or less diffuse at locations that are a priori not always easy to predict (Bour et al., 1999). 2. The stability conditions of the mining works with an outlet at the surface, such äs the shafts and adits, may be significantly and unfavourably modified in the presence of water. The treatment of these openings is generally necessary, through the implementation of a self-bearing concrete plug at the head or within the structures for example [START_REF] Wojtkowiak | Principles for a safe closure of old mine shafts and adits[END_REF], to guarantee the safety of people and property in the direct vicinity of the mine. 3. Land movements may appear at the surface consecutively to the instabilities of the bottom works and damage the frame and infrastructures. 4. Mine gas releases may also take place at the surface and the risks linked to the accumulation of these gases in confined areas and premises that are not or poorly ventilated are far from being negligible.

Based on the extent and magnitude of the geotechnical and environmental impacts on the surface of the flooding of the Underground mine works, in this articie we will only analyse the two last above-mentioned points, i.e. impacts of water rising on the stability of the surface and mine gas releases.

Firstly, we will deal with the various mechanisms that can result from the flooding of the mines and their potential effects on the surface. Based on the studies conducted by INERIS on behalf of Charbonnages de France in the various French coal basins, we will study the nature of these impacts and, if applicable, its importance vis-ä-vis the safety of people and surface facilities. Then, we will briefly present the compensatory measures the most conunonly implemented to reduce risks linked to mine gas emissions to the surface. One of these measures consists in exploiting the mine gas resource and therefore reducing the methane emissions to the atmosphere. The optimisation and the search for profitability of such an Operation, while guaranteeing the safety of people and surface facilities and the protection of the environment, lead us to develop a mathematical model that is presented in this document and that makes it possible to forecast the behaviour of the gas reservoir and ensure the perenniality of its exploitation.

GEOTECHNICAL IMPACTS OF THE WATER RISING

After Underground mines are closed, once the shafts and adits with an outlet at the surface are treated to guarantee the safety of people and property in their direct vicinity, problems linked to the long-term stability of the undermined land arise.

The rising of the Underground water, resulting from the Interruption of the dewatering pumping is one of the factors liable to modify the long-term stability conditions of the mining works, of the overburden disturbed by the works and, therefore of the surface of the land. Indeed, a change in the hydraulic environment may have a major infmence on the state of the stress forces within the massif as well äs on the mechanical characteristics of the rock mass.

Studies conducted in various coal basins made it possible to define the type of impact of the water rising on the stability of the surface and, if applicable, to translate its importance in terms of damage to the buildings.

2.l Hydraulic consequences of Underground works

While the mine is still in Operation, fractures that already exist within the massif slightly grow vertically from the works and add on to those induced by the mining operations.

All these fractures play a draining role that grants permeability coefficients to the influence volume of the mining works higher than those of the sectors located outside the area disturbed by the mine. Water from the overlying aquiferous formations may therefore flow in the direction of the mining works thus often requiring the implementation of a dewatering pumping System while the mine is in activity. The pumping gradually drives back the aquifers influenced by the activity of the mine.A desaturation of the massif complements this effect in the direct environment of the mine (dose field), whether äs regards the rock matrix or the network of fractures it is subject to.

2 Modifications of the mechanical properties of the rock mass

The state of the geological medium may be characterised by the values of a certain number of state variables pertaining to the various categories of physical phenomena, at a given moment. These phenomena are hydraulic, mechanical, thermal and chemical. This study is limited to the two first aspects (Fig. l).

The role of water on the mechanical behaviour of a massif is mainly limited to the action of the volume forces (Archimede's principle and outflow forces) and of the interstitial pressure. Each one of these forces is linked to the hydraulic load that is a function of the flooding level of the massif (piezometric altitude of the aquifers).From a mechanical standpoint, the action of water on a massif äs several consequences. • A redistribution of the stress within the massif: the transitory phase of the water rising, that is often the phase with the most unfavourable conditions for the stability of the massif, must be distinguished from the füll flooding phase that generally leads to the instalment of a new hydraulic and mechanical equilibrium within the massif. • A reduction of the mechanical characteristics of the rock mass: although it is systematic, it varies in amplitude according to the type of materials. The more the water Saturation degree is low before the flooding, the more the reduction of the geomechanical characteristics of the massif will be high. For example, we note that the ratio between the umaxial compressive strength of a dry material and that of the same material saturated in water varies from l. 2 to 2. 5 for the carboniferous sandstones and shales, that it is more or less equal to 2 for the iron ore from the Lorraine area, and may reach 4 for certain chalks from the Norm of France.

3 Consequences of the water rise on the surface stability

The Interruption of the dewatering pumping when mines are closed will gradually give rise to a rising of the Underground water that will flll residual voids that result from the mining activity until a new hydraulic and hydrogeological equilibrium is achieved.

In this case, we must plan whether or not the stability conditions of the former mining works and of the overburden that existed before the rise of the water might be significantly modified and have repercussions on the stability of the surface (subsidence restoration.local or generalised collapsing...) during the flooding. a) Case of infrastructure mining works For the mining shafts that were backfilled, the transitory flooding phase may give rise to phenomena with effects at the surface under certain conditions (untamping of the rocks, breakage of the shaft walling, breakage of the surrounding rock...). By taking into account the influence of numerous factors in favour of the stability of the rocks such äs their high granulometry or the "silo effect", it however appears that the probability of occurrence of such phenomena is low. Anyway, after the stabilisation of the water level, the risks of impact on the surface linked to the flooding become negligible or even nil for these works. Considering the geological and mining contexts of the French coal basins, the noxious influence of the rising of the Underground water on the stability of the adits can only be feit in the event of adits dose to the surface (less than twenty metres). Under these conditions, the influence of water may therefore be embodied by:

• a reduction of the mechanical properties of the overburden and of the adit supporting elements, • a reduction of the slope angle of the caving material cones consecutively to local instabilities of the roof of the adits whether reactivated or initiated by the water rising, • the influence of the hydraulic pressure.

For adits that are dose to the surface, the major risk is linked to the degradation of the mechanical properties of their supporting elements in the presence of water, which may lead to a crown-type breakage, or a further collapsing of the roof. However, after the flooding, the alteration phenomena, that may lead to a reduction of their mechanical characteristics, are considerably slowed down due to the absence of oxygen required to any alteration process.

b) The total mining works

Amongst the various types of total mining works (backfilled, caved in or drawn off longwall face), certain are more liable man others to give rise to surface effects when the Underground water rises: this is the case of the backfilled sloped longwall faces. An untamping of the rocks may take place if they are subject to three simultaneous conditions (Fig. 2). In this specific configuration, the empty volume created by the untamping of the rocks may give rise to subsidence elevations liable to open up at the surface (Fig. 3) or a restoration of surface subsidence (Fig. 4).

Forecasts conceming the magnitude of the effects on the surface of such mechanisms were performed thanks to the SUBSID Software, developed by the INERIS [START_REF] Aissaoui | Amelioration de la prevision des qffäissements mimers ä l'aide des approches empiriques[END_REF], applied to the former backfilled sloped works of the North and Pas-de-Calais coal basin. The caiculations were conducted based on the very pessimistic assumption of a systematic untamping of the rocks consecutively to the water rising. Therefore, we demonstrated that the possible subsidence restoration is negligible and may at the most slightly damage the very long buildings (over one hundred metres long all in one piece), and only if the mining works are located at a depth inferior to 50 m.

Considering the configuration of the backfilled sloped works in the basin of the North and Pas-de-Calais and the type and thickness of the overburden, we also highiighted the fact that the risk of subsidence elevation should be excluded.

c) Partial mining works

The major instability risks are linked to the partial mining works referred to äs mining works by rooms and abandoned pillars (Vinkler & Piguet,1999). The alteration of the mechanical characteristics and gradual distortions that affect the various layers, in particular those of the roof, and the pillars have sometimes resulted in the failure of former works several dozen years after they were abandoned.

To this regard, the action of the water is unfavourable to their stability in the short term, äs it may give rise to a collapse or subsidence. On the contrary, in the long term, once the layers have retumed to their initial piezometric level, no Underground extemal factor can influence the stability of the massif. However, when the land is sensitive to creep that the flooding can reactivate, the Situation may change over several years until it is stabilised or until instability is triggered.

4 Elevation of the topographical surface

In certain foreign coal basins, the rising of the undergroimd water is accompanied by an elevation of the surface in the area influenced by the former Underground nüning works (Bekendam & Pöttgens, 1995). The mechanisms at the origin of this phenomenon still form the subject of studies, however, they seem to be linked to a reduction of the effective stress. In any case, the amplitude of the elevation measured at a given point is much lower than the amplitude of the subsidence supported by the same point during the mining activity. The elevation only represents several percent of the amplitude of the total subsidence. Furfhermore, this elevation is gradual and concems larger areas comparable to those of the subsidence basins. Therefore, there is not much chance that this type of elevation will damage the frame and infrastructures on the surface. Indeed, damage is essentially linked to the horizontal deformation of the ground that is very low in this case.

The levelling measurements performed in several French coal basins that are currently being flooded, will constitute a knowledge and database unique in France to further detail the phenomenon in terms of amplitude and duration, äs well äs the mechanism(s) at the origin of this elevation.

Conclusions

In most of the studied total mining works, the effects on the surface linked to the water rising are in fact only a restoration of land movements due to the reactivation of mechanisms already initiated by and during the Underground mining activity. Most of the time, their amplitude and duration are extremely moderate. Furthermore, the fact that most of the risks of impact on the stability of the surface are much higher during the transitory phase of the flooding, should be highiighted. New mechanical balance conditions are set within the undermined massif with a retum to a hydraulic and hydrogeological balance of the aquifers. Finally, knowledge äs regards the role that the interstitial pressure plays on the long term behaviour of the rocks is still very poor. This partly explains why the hydromechanical coupling is only slightly or even not at all taken into account in the numerical modelling of the behaviour of the Underground works in fractured rock massifs.

RISKS AND PRECAUTIONARY MEASURES LINKED TO GASES IN CLOSED MINES

3.l

Firedamp accumulation mechanisms and gas emission to the surface After the closure of a coal mine, firedamp emissions can stop in certain cases but some gas remains in the void spaces of the Underground mine (mining works, fractures and cracks of the rock mass...). Gas emissions can also continue for a while after the closure. However, the gas emissions can stop with the flooding of coal layers, up to several meters above the coal seams. Therefore, the gas accumulated could migrate to the surface if the pressure within the reservoir is superior to the atmospheric pressure at the surface.

This difference of pressure between the reservoir and the surface may be caused by:

1. an increase of the pressure inside the reservoir with :

• the water rising within the Underground void spaces that compress the gas,

• the gas firedamp desorption from the coal in place; 2. the decrease of the atmospheric pressure due to meteorological depressions. Furthermore, all these possible causes of gas emissions can take place at the same time. Gas then can migrate to the surface : 1. in a privileged manner through communication works between the mine and the surface (old shafts, adits...) ; 2. through the overburden itself.

2. Type of risk linked to the firedamp

The firedamp in the mining caviting and migrating to the surface can cause two risks: • an inflammation risk, because methane constitutes a flanunable mixture with air in certain proportions (from 5 to 15 %): • a risk of asphyxia, because firedamp accumulates in places not or badly ventilated replacing the atmospheric air. Several examples of traditional hazardous situations at the origin of incidents or accidents are indicated below: -punctual and concentrated firedamp emission through a direct communication between the former works and the surface (mining opening, crack...); -accumulation of gas in confined or semi-confmed areas (cellars or basements, even houses or buildings themselves, but also Underground networks...); -restoration of the communication of old works with the atmosphere during drillings or civil engineering works....

3 Prevention and monitoring measures

Beyond the necessity to treat the mining openings, several compensatory measures (treatments) may be implemented to limit the migration of the fire-damp and its aecumulation up to explosive contents. Two types of measures may be defined: -the passive treatment: in certain cases, it may be necessary to perform bore holes from the surface to intercept certain high points of the Underground reservoir. On the one hand, they make it possible to check the parameters of the gas reservoir : composition and pressure. On the other hand, they can constitute outlets for the gas that can be totally controlled and safe (far from housing areas and fitted with safety devices such äs shut-off valves), and can more or less partially replace the noncontrolled "natural" outlets.

-the active treatment: the Underground reservoir can be in depression with the atmosphere by extracting the gas. If surface communications are sufficiently tight, it is therefore possible to reach very low absolute pressures from 0. 3 to 0. 5 bar. The extent of influence of a drainage borehole can be very large: several km2 and even hundreds of km2.

4. Drainage installations

In French coal basins, gas drainage installations on abandoned mines appeared in response to the following objectives.

• Safety: the mine gas drainage can be an efficient way of controling gas emissions to the surface. • A contribution to the protection of the environment, by the reduction of the methane emissions to the atmosphere. • Energy recovery, with the marketing of the gas.

Facing the requirements to extract and market such a resource, modelling tools are needed to forecast or evaluate the behaviour of these gas reservoirs. As they did not exist, a mathemaäcal model was therefore developed to appreciate the capability of the reservoir to supply gas for a given period, and to characterise the operating conditions that will guarantee an Optimum production. a)

Characteristics of a mine gas reservoir b)

A firedamp reservoir is the entire volume destabilized by an old mining (a district, a main level, or a mine), made tight from the outside and whose voids spaces grew rieh out of methane.It is often about a large-sized deposit, whose gross volume can reach several tens of million cubic meters. Contrary to a conventional reservoir, volume is not materialized by defined walls or extemal immediate strata. üsefül volume corresponds then to the volume of the cavities and spaces due to the cracking and the porosity of remaining coal äs well äs to residual volumes of the galleries and the Underground workings. This volume does not represent the mine gas reserve since the methane is at the same time stored in a free form in the pores and cracks of the coal and rock matrix, and under an adsorbed form on the coal surface within the micropore structure. In case of communication with a water table, the reservoir volume decreases progressively with the rising of water. b) General behaviour of the reservoir In extracting mode, the outgoing flow of drained gas is generally higher than the incoming flows of gas within the reservoir. Then, the reservoir pressure decreases, which increases parasitic air intake in the reservoir. The gas drainage can thus be continued until a certain absolute reservoir pressure : from a pressure level, the incoming flow of air involves in general a too important reduction in the methane content of the drained gas, which must be at least 50% to 65%. Periods during which the extracäon is stopped must thus be programmed in order to reload the reservoir with methane.The absolute pressure within the reservoir and the methane content of the drained gas are the two principal parameters goveming the behaviour of the reservoir. The figure 6 presents, for a real site, the evolution of the pressure within the reservoirs and the methane content of the mine drainage. c)

Reservoir modelling A modelling tool was developed to forecast the behaviour of the reservoir (pressure within the reservoir and methane content in the collected gas) according to the extracted gas quantities.

Due to the coal structure, that may be assimilated to an assembly of grains drained by a network of cracks, the study of mine gas storage takes place at two levels:

• äs regards the material gram where the adsorption, desorption and diffüsion phenomena take place, • äs regards the network of cracks in which the "free" gas flows out according to the permeability. Therefore, the physical laws goveming the gas circulation in the reservoir are complex:

• firedamp desorpäon (Langmuir law...), • mine gas diffüsion within a coal massif (Fick law...), • gas outflow in porous and permeable media (Darcy law...),

• circulation of fluids in niining infrastmcture works and gas pipes.

The approach consisted in comparing the reservoir to a total volume whose pressure varies with the outgoing flow of mine drainage gas and the total incoming flow due to the firedamp desorption and the air intake.The law of perfect gases was retained to connect the quantity of gas contained in the reservoir and the pressure in this one. For relatively steady conditions, it is then possible to caiculate the total incoming flow of gas (firedamp desorption plus air intake) and the volume of the reservoir.

The modelling suggested was implemented using experimental data from two reservoirs currently exploited and located m the Lorraine basin in the East of France.The model was evaluated on the results of forcasts based on past experimental data (Fig. 6). The results prove very satisfactory with, for the modelled variables (pressure and content mefhane), an absolute average error between the observed and predicted values lower than 10%. Discrepancies between observed and predicted values can also be observed during the stop of the extraction : the absence of gas circulation makes the measurements much less representative of the gas reservoir. Although the gas sources localization remains difficult, one can indeed easily admit that gas drainage, air intake and desorption can take place in three distinct points in the reservoir. The measurements of methane content is made at the level surface, on the gas extraction installation. Due to the large size of the deposit, the measurements representativeness remains relative.

Results

The pressure of the reservoir and the methane content of the extracted gas have been modelled according to the quantities of extracted gas. The periodic caiculation of the model parameters with recent exploitation data will allow a follow-up of the site characteristics in terms of reservoir volume and desorption potential.

Lastly, modelling must also be employed in order to optimize the phases between the stops and the extraction, or the mine drainage gas flowrate.

The model thus developed and validated on two real cases could be used thereafter on other sites and will constitute a true tool of optimization and industrial management of a firedamp reservoir.

CONCLUSION

Based on numerous studies performed in the various French coal basins, the INERIS developed a comprehensive method to identify and analyse the risks linked to the definitive closure of imderground mining works, especially by integrating the consequences on the surface of the water rising within the former mining works. Grounded on this analysis, compensatory measures to reduce these risks may be defined and optimised whether from the technical or economic standpoint, while guaranteeing the safety of people and property in the long term. Over time, the implementation of measurement and monitoring devices in the closed mining basins should provide data indispensable to compare the forecast with the actual behaviour of the mining works.

ACKNOWLEDGEMENTS

The authors would like to thanks Charbonnnages de France for their technical and financial support to perform the risk studies in the French coal basins.

Legend of Figure 2: The three simultaneous conditions for the backfill outflow of a stowed sloped face.

Legend of Figure 3: Sinkhole mechanism initiated by the backfill outflow of the stowed sloped face with the presence of Wealdian sands in the overburden (North and Pas-de-Calais coal basin).

Legend of


Frederic Y Bois
email: frederic.bois@ineris.fr

Dr Frédéric Y Bois

Applications of Population Approaches in Toxicology

Keywords: Bayes, Toxicodynamic models, Toxicokinetic models, Population Models

Many experimental or observational studies in toxicology are best analysed in a population framework.

Recent examples include investigations of the extent and origin of intra-individual variability in toxicity studies, incorporation of genotypic information to address intra-individual variability, optimal design of experiments, and extension of toxicokinetic modelling to the analysis of biomarker studies in industrial hygiene or epidemiological contexts. Bayesian statistics provide powerful numerical methods for fitting population models, particularly when complex mechanistic toxicokinetic or toxicodynamic models are involved. Challenges and limitations to the use of population models, in terms of basic structure, computational burden, ease of implementation and data accessibility, are identified and discussed.

INTRODUCTION

Uncertainty and variability affect potentially all experimental results. Both are prevalent when dealing with biological material or with human subjects. While uncertainty can be construed as a nuisance to Cartesian minds, variability has a different status. It may be necessary to assess both, for example when predicting risks. In such case, confidence in predictions and evaluation of the range of individual risks within a population are important considerations. Uncertainty is clearly different from variability:

-Uncertainty is essentially due to lack of knowledge [START_REF] Rowe | Understanding uncertainty[END_REF] and may have various sources.

Measurements are made only with finite precision. When data are used to fit toxicokinetic (TK) or toxicodynamic (TD) models, there is always some "noise" about the resulting model parameters values. In addition, the fact that a population sample is studied also introduces an element of randomness. Simplifications or mispecifications of the fitted models also yield parameter uncertainties. Finally, extrapolations of in vitro results to in vivo, or across species, etc., are also a large source of variance.

-Variability refers typically to differences in characteristics among individuals (inter-individual variability), or across time within a given individual (inter-occasion, or intra-individual variability).

It may stem from genetic differences, lifestyles, physiological status, age, etc. [START_REF] Rowland | Clinical Pharmacokinetics: Concepts and Applications[END_REF] . Experimental toxicology studies are in general specific to the (few) individuals or animal studied and extrapolating their results to larger and different populations requires careful consideration. TK studies for example are rarely about the toxicokinetics of a compound in a particular individual (clinical pharmacokinetics may be different in that respect). They rather aim at information about average toxicokinetics in sensitive populations, age groups, whole country population, etc., and at a quantification of the variability in those kinetics.

Uncertainty may be reduced by new experiments, the design of which can be optimised to bring a maximum of information [START_REF] Beatty | Optimal statistical design for toxicokinetic studies[END_REF] , or by a better understanding of the biology and hence better models. This may take much efforts and decisions routinely need to be quickly made despite a fair amount of uncertainty about the situation examined. The "precautionary principle" is indeed an answer to the recognition of that need [START_REF] Stirling | On science and precaution in the management of technological risk -Final report and summary of the study conducted under the auspices of an ESTO project on technological risk and the management of uncertainty[END_REF] . Variability is inherent in animal and human populations. It cannot be reduced and seems even desirable for the stability of the global ecosystem [START_REF] Mccann | The diversity-stability debate[END_REF] . The fact that uncertainty is almost always present and hampers a precise assessment of variability, has led to the development of somewhat sophisticated statistical "population" models to disentangle them. A brief account of the development of those tools, of their basics, and a short survey of their burgeoning applications to toxicology form the topic of next sections.

POPULATION MODELS

Population models were originally introduced in the context of pharmacokinetic studies for drug development and evaluation [START_REF] Beal | The NONMEM system[END_REF][START_REF] Beal | Estimating population kinetics[END_REF][START_REF] Sheiner | The population approach to pharmacokinetic data analysis: rationale and standard data analysis methods[END_REF] . A number of applications can be found in the pharmacokinetic literature [START_REF] Aarons | Population pharmacokinetics[END_REF][START_REF] Sheiner | Population pharmacokinetics/dynamics[END_REF][START_REF] Vozeh | The use of population pharmacokinetics in drug development[END_REF][START_REF] Yuh | Population pharmacokinetic/pharmacodynamic methodology and applications: a bibliography[END_REF] . Until recently, little attention had been paid to these approaches in toxicokinetics. This is partly due to the fact that toxicokinetics used a lesser amount of human data, and focused on animal models. That picture is changing with the rapid development of human biomarker studies and improvements in analytical chemistry, which make is possible to follow in humans the kinetics of trace amounts of chemicals. Before our work in this area (Bois et al., 1996a;Gelman et al., 1996) , the only published attempt to link population variability to TK modelling was presented by [START_REF] Droz | Variability in biological monitoring of organic solvent exposure -II. Application of a population physiological model[END_REF] , but no fitting to experimental data was performed in that exercise. Indeed, the scope of population models is not limited to pharmaco-or toxicokinetics. Dose-response studies, for example, can be fruitfully analysed with the same tools [START_REF] Davidian | Nonlinear Models for Repeated Measurements Data[END_REF][START_REF] Palmer | Bayesian optimal design in population models for haematologic data[END_REF][START_REF] Patron-Bizet | Assessment of the global two-stage method to EC 50 determination[END_REF][START_REF] Racine-Poon | A Bayesian approach to nonlinear random effects models[END_REF][START_REF] Racine-Poon | Population models[END_REF][START_REF] Wakefield | Bayesian analysis of linear and nonlinear population models by using the Gibbs sampler[END_REF] .

The objective of population models is to extract, from data on individuals units (e.g., subject, family, group), a quantitative description of the variability of the kinetics of a compound within a large population.

The basic idea is that the same model can describe the data for each unit, and that the model parameters differ randomly between units, or between occasions within the same unit, etc. Such randomness characterises variability per se, and is described by (multivariate) probability distributions. Individuals' renal clearance, for example, will be assumed to be normally distributed around a "population mean", with a "population variance" which measures variability in the population. Population means and variances (one for each kinetic parameter supposed to vary across units) are termed "population" parameters. They are estimated by data fitting, together with the model parameters for each unit. Other population parameters could be defined, such as parameter covariances (if some parameters are supposed to be correlated at the population level), or degrees of freedom if a (student) t-distribution is assumed instead of the normal distribution.

A typical population TK model is presented in Figure 1 (Bois et al., 1996a;Bois et al., 1996b;Gelman et al., 1996;[START_REF] Wakefield | The Bayesian analysis of population pharmacokinetic models[END_REF] . The "hierarchical" or "multilevel" structure of the model is apparent, with two major components: the subject level and the population level.

At the subject level, data (y) are measured at times t, after exposure to dose E. A set of individual characteristics (parameters, ) also condition the observations. Some of these characteristics maybe known (e.g., sex, age, genotype), in which case they are considered as "covariates", but usually some interesting characteristics are unknown and need to be estimated from the data. That estimation can be achieved by data fitting. In Figure 1, unknown quantities are enclosed in circles and known quantities in square boxes. Data fitting requires basically the definition of two models:

-First, a statistical model should link some "predictable" aspects of the data (e.g., their expected values, or their variance) to the actual measurements (y). Deviations between measurements and their expected values may be due to limited assay precision, modelling error, or random intraindividual variability. Such deviations are typically assumed to be independent normal or lognormal random variables, with mean zero and variance  2 . That variance can be set to a known value (e.g., if analytical precision has been independently assessed and if modelling error is negligible). It can also be estimated along with the other parameters. The model chosen will specify the data "likelihood". That likelihood is proportional to the probability of the observed data, given a model and its parameter values.

-Second, a TK or TD model, f, should link the chosen "predictable" to E, t, and , or any other relevant "explanatory" variable.

At the population level, inter-unit variability is described by considering that the individual parameters, , are randomly drawn from given statistical distributions, each with population mean µ and population variance ∑ 2 , as described above. Population means and variances are to be estimated together with the parameters  and the residual variance  2 . What about uncertainty in the above framework? Uncertainty enters in the fitting process in two ways:

-First, through all variance sources aggregated into the residual variance  2 (typically, limited assay precision, modelling error, and un-modelled levels of variability). Note that it is possible to disaggregate  2 in order to identify separately its various components. For example, a further level in the hierarchy can be defined to deal with intra-unit variability [START_REF] Lunn | Markov chain Monte Carlo techniques for studying interoccasion and intersubject variability: application to pharmacokinetic data[END_REF] . In that case, residual variance can be reduced if intra-unit variability is significant. However, some measurement and modelling errors are likely to remain present, even with the best available data and models. A side-effect of that "noise" is to allow different model parameter values to yield adequate fits to the data (what is meant by adequate is determined by the statistical measurement model).

That results in uncertainty about all the quantities estimated in the fitting process (population parameters, individual characteristics, etc., and even  2 if it is adjusted). This uncertainty may entangle estimates through fitting-induced correlations, which make it difficult to interpret the results when many parameters are adjusted.

-A less obvious source of uncertainty, which affects the population parameters, is the finite and often small size of the population sample studied. Even if the characteristics of each individual were perfectly known (from perfect data and models), the population parameter values inferred from a non-exhaustive sample would still be affected by uncertainty.

What are the advantages of population modelling? Most of all, it is a clean way to model levels of variability in data obtained on populations. The "naive" approach which consists in fitting individual data separately, collecting individual parameter estimates and forming their average etc., does not work if very precise parameter estimates cannot be obtained. In clinical pharmacokinetics and in toxicokinetics, the data on each individuals tend to be sparse, for ethical, feasibility or cost reasons. Such data usually lead to fairly uncertain parameter estimates and population approaches are then warranted. In addition, in a population model, information of each subject is reinforced by "borrowing strength" from the other subjects' data and the estimation process is improved overall. For example, when a subject's characteristic is barely estimable with its own data, the TK or TD link model can fail to give any reasonable estimate for the other characteristics (an effect of fitting-induced correlations, as in the flip-flop effect in pharmacokinetics).

Discarding subjects data for that reason can make the experiment or trial worthless. In a population model, each subject's characteristics will tend to resemble that of the population (subjects are drawn toward the population mean), in particular if little data is available on that subject, and the estimation can proceed without information loss.

BAYESIAN ANALYSIS OF POPULATION MODELS

A number of methods have been proposed to fit models and estimate their parameters within a population framework. One way to classify them is into parametric and nonparametric methods. In the parametric approach, the shape of the population distribution is assumed to be known, and the population parameter values are quantities to estimate [START_REF] Beal | Estimating population kinetics[END_REF][START_REF] Racine-Poon | Population models[END_REF] . That approach is the most often used. In the nonparametric case, no assumptions are made about the shape of the population distributions, and both shape and parameters are estimated from the data [START_REF] Davidian | Smooth nonparametric maximum likelihood estimation for population pharmacokinetics, with application to quinidine[END_REF][START_REF] Mallet | Nonparametric maximum likelihood estimation for population pharmacokinetics, with application to cyclosporine[END_REF] .

Another classification of the relevant statistical methods differentiates maximum likelihood [START_REF] Beal | The NONMEM system[END_REF] and Bayesian approaches [START_REF] Racine-Poon | A Bayesian approach to nonlinear random effects models[END_REF] . Modelling in toxicology ought to be particularly well served by a Bayesian approach. TK/TD models tend to deal with scarcity of data and extrapolation needs (e.g., inter-species, inter-dose) by relying on complex mechanistic structures [START_REF] Kohn | Achieving credibility in risk assessment models[END_REF] . Physiologically-based toxicokinetic (PBTK) models [START_REF] Bois | Analysis of PBPK models for risk characterization[END_REF] or cancer models [START_REF] Tan | Stochastic Models of Carcinogenesis[END_REF] can be fairly sophisticated, compared to regression or classical compartmental models. In a Bayesian framework, nonlinear models and non-standard statistical distributions are easily dealt with [START_REF] Bernardo | Bayesian Theory[END_REF] . More importantly, a priori physiological information (such as found in the literature, in vitro, etc.) about parameters can be seamlessly incorporated into the data analysis (Gelman et al., 1996) . Finally, modern numerical techniques have made Bayesians unafraid of the many parameters of population models.

In a Bayesian setting, all unknowns, such as , µ and ∑ 2 in the above example (Figure 1), are considered as random variables. The probability distribution of an unknown parameter is interpreted in terms of degrees of belief about its possible values. Before analysing an experiment, a "prior" probability distribution, P, is constructed for each unknown to reflect current knowledge about it. The use of a probability distribution to represent prior knowledge offers a unified way to account for precise as well as vague information available before an experiment. When little knowledge is available on a particular parameter, a non-informative distribution (i.e., flat-shaped) can be used. If stronger prior knowledge is at hand, the prior should represent that information as appropriately as possible. Note that information from the literature tends to relate to average values or variances for a population and is naturally handled in a population model. With a Bayesian approach, such information can be directly used under the form of prior distributions for population parameters such as µ and ∑ 2 (see above).

The set of prior distributions is then updated using the data values, y, to yield a "posterior" probability distribution. While the prior distributions may be specified independently from each other, the posterior is a joint distribution of all unknowns (i.e., correlations induced by the model fitting process are accounted for).

Bayes' theorem indicates that the posterior probability distribution of , µ and ∑ 2 given y, noted P(,µ,∑ 2 |y), is proportional to the product P of all priors by the data likelihood defined above, P(y|,µ,∑ 2 ):

P(,µ,∑ 2 |y)  P  P(y|,µ,∑ 2 ) (1)
The hierarchical structure of a population model allows a splitting of the likelihood in two simpler terms:

P(,µ,∑ 2 |y)  P  P(y|) P(|µ,∑ 2 ) (2)
The exact form of term P(y|) is dictated by the error model and the TK/TD model chosen. The form of P(|µ,∑ 2 ) depends on the distribution assumed for the individuals' characteristics within the population.

Bayes' rule can be viewed as a formula that shows how existing beliefs, formally expressed as probability distributions, should be modified by new information [START_REF] Bernardo | Bayesian Theory[END_REF] . Any statistical summary (mean, standard deviation, quantiles, correlations, marginal distribution, etc.) about the model parameters can be obtained from their posterior distribution. However, it is often difficult to obtain an analytical expression for P(,µ,∑ 2 |y). Until recently, that has been a strong impediment to practical Bayesian analyses.

Fortunately, the advent of Markov chain Monte Carlo (MCMC) methods has alleviated a number of these difficulties. These methods are powerful tools which can provide samples of parameter values out of their posterior distribution even without knowledge of an analytical expression for it [START_REF] Gelman | Markov chain Monte Carlo methods in biostatistics[END_REF][START_REF] Gilks | Markov Chain Monte Carlo in Practice[END_REF] . They have largely contributed to the recent proliferation of Bayesian analyses in applied statistics.

APPLICATIONS

TOXICOKINETICS

Since population concepts were ported to toxicology via toxicokinetics, most applications seen so far fall into this section. TK models can be used to relate external exposures to toxicants to internal measures of dose. They typically describe the body as a set of compartments through which chemicals travel or are transformed. TK models fall into two main classes: "classical" models and physiologically-based (PBTK) models.

Classical TK models represent the body by several connected compartments (rarely more than three) within which chemicals are assumed to distribute homogeneously [START_REF] Gibaldi | Pharmacokinetics[END_REF] . Compartments do not necessarily reflect the anatomy of the species of interest. The number of compartments and parameter values (e.g., exchange rates between compartments, volumes) are inferred from fitting the model to TK data. These models are thus commonly called empirical, or data-based. Population models can be used which such models as link functions between the data and individual characteristics, as is routinely done in pharmacokinetics.

In PBTK models, the body is subdivided into a series of anatomical or physiological compartments [START_REF] Bailer | An introduction to the use of physiologically based pharmacokinetic models in risk assessment[END_REF] . Connections between compartments represent the blood or lymphatic circulation. PBTK models typically use three types of parameters: physiological parameters, such as breathing rate, blood flows and tissue volumes; physicochemical parameters, such as partition coefficients; and biochemical parameters, describing for example metabolic processes. The number of compartments depends on the objective of the study, and of the possible mode(s) of action and site(s) of toxicity of the chemical studied. The focus of PBTK modelling has been mainly on improving the scientific basis of animals to humans extrapolations [START_REF] Gargas | Pharmacokinetic modeling approaches for describing the uptake, systemic distribution, and disposition of inhaled chemicals[END_REF] . The ability to tailor the parameter values of these models to a particular animal species offers, in theory, a firm basis for interspecies extrapolation. In fact, for most PBTK models, some parameters (in particular those describing metabolism) are difficult to extrapolate across species. A partial answers to this question may be provided by the "parallelogram approach" [START_REF] Morris | A parallelogram approach for safety evaluation of ingested acetaldehyde[END_REF] , but more often, some parameters values need to be obtained through model fitting to TK data. Obviously, that fitting should take into account the available physiological prior information on parameters. Prior information, as mentioned above, is usually about population means or variances, and population modelling offers a natural setting for inclusion of such information. To further complicate matters, PBTK models have many parameters (10 to 20 is typical) and are often nonlinear and without analytical solution. It is only recently that a general Bayesian method has been presented for fitting these models (Bois et al., 1996a;Gelman et al., 1996) . Before then, it was impossible to decide whether the fits were acceptable, the models reasonable, and what confidence to place in their results. A treatment of these models in a Bayesian population framework can therefore be advantageous, as further illustrated below, where we focus on specific study types.

Applications to Experimental Studies

Animal Studies

The analysis of animal experiments may benefit from the use of population models. Inter-animal variability in kinetics exists as well as in humans. It may not be so high because of the use of inbred strains, but is still there. Other levels of variability can also be present, such as inter-lot or inter-strain variability. If a large number of data points is gathered on each experimental unit and if characteristics of interest can be precisely assessed, a naive approach might be sufficient. However, when a PBTK model is used, the precision of estimates obtained is not guaranteed, given the number of parameters in those models. Using "standard" values for most physiological parameters and adjusting just a few "sweeps uncertainty under the rug" and leads to over-optimistic estimates of precision. A careful analysis would therefore proceed along the population modelling way, unless it can be proven that a simpler approach is also correct.

Human Volunteers Studies

Studies in human volunteers are being performed, despite potential ethical concerns, for several reasons; i.

Modern analytical techniques allow the determination of kinetics after very low exposures (within environmental or occupational exposure ranges) posing minimal or zero risk to the subjects; ii. Animal studies may not be predictive of human responses; iii. Ethical concerns also exist for animal studies, even if rodents do not get to sit on ethics committees. Human toxicology studies are clearly geared toward risk assessment rather than toward research on mechanism of action, for example. In that context, an independent assessment of variability is a clear concern, and population modelling is a natural answer.

PBTK models are here also in favour, since they can be used to simulate toxicokinetics under special conditions, such as physical activity [START_REF] Dankovic | The impact of exercise and intersubject variability on dose estimates for dichloromethane derived from a physiologically based pharmacokinetic model[END_REF] , or in specific targeted populations (e.g., pregnant women, children), by simply changing the values of some physiologic parameters. The same arguments as for animal experiments can be made in favour or the use of population models in that context.

For further discussion, the reader is referred to our recent population TK analyses of volunteer exposures to tetrachloroethylene (Bois et al., 1996a;Gelman et al., 1996) , benzene (Bois et al., 1996b) , trichloroethylene (Bois, 2000b;Bois, 2000a) , or 1,3-butadiene [START_REF] Bois | Optimal design for a study of butadiene toxicokinetics in humans[END_REF] . The inclusion of genotype information in such analyses is straightforward [START_REF] Johanson | Development of new technique for risk assessment using physiologically based toxicokinetic models[END_REF] , and should see considerable development in the near future.

Optimal Experimental Design

A useful extension of data analysis is the optimal design of experiments. Population modelling and Bayesian statistics can be jointly used to determine, for example, the optimal sampling times and optimal number of subjects to study [START_REF] Atkinson | Optimum experimental designs for properties of a compartmental model[END_REF][START_REF] Merlé | Bayesian design criteria: computation, comparison, and application to a pharmacokinetic and a pharmacodynamic model[END_REF][START_REF] Wakefield | An expected loss approach to the design of dosage regimens via sampling-based methods[END_REF] . In toxicokinetics, we recently used a classical two-compartment model to address similar questions for the planning of human volunteer exposures to 1,3-butadiene [START_REF] Bois | Optimal design for a study of butadiene toxicokinetics in humans[END_REF] . We also recently established that numerical sampling techniques are powerful enough to allow the use of a PBTK model in the context of design optimisation (unpublished results).

Applications to Biomarker Studies

Biomarker studies of exposure (or effects) might well be the best field of application for population TK modelling. Multivariate endpoints are often measured and it makes sense to model their temporal correlations or autocorrelations via interpretable and knowledge-based TK models rather than through empirical covariance models. The small number of measurements usually made on individuals and the medium to large size populations investigated make them potential case-studies for the population approach. However, at least three reasons have slowed development in that area:

-First and foremost, biomarker studies (unless they use controlled exposure conditions, in which case they would fall under the "experiment" heading) have uncontrolled exposure conditions. It is possible, in a population statistical framework, to consider exposure as one of the parameters to estimate. This estimation challenge is known as the error-in-variables problem [START_REF] Seber | Nonlinear Regression[END_REF] . Exposure reconstruction is a step in that direction [START_REF] Georgopoulos | Reconstruction of short-term multi-route exposure to volatile organic compounds using physiologically based pharmacokinetic models[END_REF] , although it is not always well statistically founded. A major problem, however, resides in accounting fully for the uncertainties stemming from unknown time-varying exposures. The impact of particular functional form for the time evolution of exposure has not yet been thoroughly studied and validated.

-Second, it is difficult to find a balance between "nice", powerful, physiological, complex PBPK models and simple, economical, easy to fit classical models when the number of individual studied gets large. This may well be a purely computational problem, but we are facing it. Minimal (or reduced) PBTK models may be promising tools to analyse data where many subjects are involved.

In these models, some of the compartments (e.g., poorly and well-perfused tissues) have been lumped together. Reducing the number of compartments can be a way to simplify the model while partly preserving physiological realism. In the case of benzene, we showed that a reduced threecompartment semi-PBTK model can be a good alternative to a standard five-compartment PBTK model [START_REF] Woodruff | Structure and parameterization of pharmacokinetic models: their impact on model predictions[END_REF] . [START_REF] Thomaseth | Estimation of occupational exposure to 2,3,7,8-tetrachlorodibenzo-pdioxin using a minimal physiologic toxicokinetic model[END_REF] considered a minimal physiological model to describe long-term kinetics of tetrachlorodibenzo-p-dioxin in human tissues, using data from a cohort of 359 subjects. This approach also needs to be further validated, but we feel that, as any modelling exercise, the development of a TK model should be essentially guided by the question raised and the data available to answer it.

-Third, the unfamiliarity with the approach. Here the situation is definitely improving and, hopefully, articles like this one will contribute to a better appreciation of the approach.

TOXICODYNAMICS

All the above remarks apply equally well to TD modelling, even though populations approaches are less developed in that area than in TK modelling (mostly for historical reasons). Biomarkers of effects, for example, call for the same treatment as biomarkers of exposure. Dose-response models can easy replace TK models in the population framework. An interesting application, at the TK/TD interface, is found at the occasion of drug exposure assessment in preclinical toxicological studies. In small rodents, repeated sampling may alter the health of the animals and jeopardise toxicity evaluation. To circumvent this limitation it is possible to collect serial samples from satellite animals maintained in conditions similar to those of the main study animals but not monitored for toxicity. Main study animals can then be only sparsely sampled. As noted above, sparse sampling designs are best analysed in population framework. [START_REF] Burtin | Sparse sampling for assessment of drug exposure in toxicological studies[END_REF] evaluated this option and showed that the data, analysed with a classical TK population model, yielded adequate estimation of drug disposition in both groups. An interesting extension would be to combine PBTK and PD models to analyse biomarker/epidemiological data on the toxicity of chemicals.

PERSPECTIVES

This rapid survey of the applications of population modelling to toxicology would not be complete without a few words about its limitations. The foremost may well be about a basic assumption of these models. To what extent is it correct to assign to randomness what could be, in fact, determined by time varying factors or other causal effects? When inter-individual differences could be due to age, for example, wouldn't it be better to try to model age as an explanatory variable? A very deterministic vision of the world would assume that the population approach is just a temporary band-aid while waiting for better models. It is likely, however, that population models will play an increasingly important role in exploratory analyses of heterogeneity of response to toxic exposures. They can at least point to research directions about the determinants of important differences in kinetics and actions of toxicants.

Another set of question is posed by the consideration of fully stochastic models, such as probabilistic cancer models [START_REF] Bois | Sensitivity analysis of a new model of carcinogenesis[END_REF][START_REF] Sherman | Calculation of the cumulative distribution function of the time to a small observable tumor[END_REF][START_REF] Tan | Stochastic Models of Carcinogenesis[END_REF] . Stochastic models can explain heterogeneity in response by the random evolution of a process, without appealing to the hypothetical variability of some parameters. They tend, however, to be very complex and difficult to fit to data. It would be very interesting to associate population and stochastic modelling, in the case of cancer, for example (where the respective roles of "bad luck" and genetic makeup are not easy to dissociate).

A current priority could be the development of friendly, fast, flexible and cheap software for Bayesian population modelling. It is currently difficult to quickly explore alternative TK/TD models if they include more than about 5 parameters. Computer power increases daily, but that provide only part of the solution to the problem. Life is also too hard for Bayesians, for they routinely need to define prior distributions. That is not always an easy task, even when it would seem obvious, for reasons of data accessibility. While it is well known that some TK parameters are likely to exhibit inter-or intraindividual variability in humans (or other species), the only values readily available are "reference" values for young Caucasian males. Prior information about population variability may eventually be found by going back to the original publications but even those often lack adequate statistical treatment. Data about the shape of the distributions are even harder to find. That situation really needs to be improved. Despite these limitations, experience has been sufficient to demonstrate the feasibility and usefulness of the population approach in toxicology. Its potentials are just awaiting fuller exploitation.

FIGURE CAPTIONS

Figure 1 :Figure 1 :

 11 Figure 1: Graphical representation of a population model. At the individual level, exposure (E), time (t) and individual characteristics (parameters, ) condition the observations (data, y). The relationships between E, t,  and y are described by a model, f. Individual characteristics are randomly distributed in the population. Population means and variances (µ and ∑ 2 ) describe the shape of that distribution. Measurement errors, modelling errors, and intra-individual variability are lumped in a common variance term  2 . Unknown quantities are in circles, known quantities in square boxes.


M Sylla

M Ai Heib

J P Piguet

3D numerical modelling of underground excavations in a faulted rock mass using the Boundary Elements Method (BEM) Applications to the Provence colliery

This paper presents a 3D-numerical modelling technique for underground excavations in a faulted rock mass. The displacement discontinuity method is used to solve the differential equation of the problem. A presentation of this method is given in the first part of the paper. The second aspect of this paper is concerned with the verification of the new code by comparison with the results obtained from commercial numerical codes UDEC (distincts elements) and CESAR (finites elements). Case study involving a coal mining panel excavation in a faulted rock mass is presented to demonstrate the application of the code to practical mine problems.

Résumé :

Ce papier présente le développement d'un code de modélisation numérique permettant de calculer les contraintes induites par une exploitation souterraine en présence de failles. La méthode des discontinuités de déplacement est utilisée pour résoudre l'équation différentielle du problème. Une présentation de cette méthode et son utilisation pour modéliser des veines et des failles sont présentées dans la première partie de cet article. Ensuite, une vérification de ce code, par comparaison aux codes de calculs UDEC (éléments distincts) et CEASR (élément finis), est proposée. Enfin, une étude de cas portant sur l'exploitation d'un nouveau panneau en présence de failles dans la mine de Provence (Sud de la France) est mené pour confronter les résultats obtenus à un cas réel.

Introduction :

The presence of faults in the rock mass surrounding underground mines can cause failure and even major rockbursts due to sliding on fault planes. Therefore, for safety and economic reason, it is necessary to be able to predict the behaviour of a mining zone, particularly, when it is situated in a faulted area. To achieve this goal, the authors have recently developed a 3D numerical modelling code (FAULT3D) using the Boundary Element Method (BEM).

Model input parameters include rock mass parameters (Young modulus, Poisson's ratio, unit weight, initial or in situ stresses), fault characteristics (tensile strength, cohesion, friction angle, normal and shear stiffness), mine geometry (panels and faults). A method to treat underground voids (gob or backfill) was developed. A 2 or 3 dimensional geometry view is available.

The BEM requires to calculate and to solve a system of equations resulting from the influence of all the boundary elements. That is quite a big advantage compared to an equivalent differential method but it rapidly becomes inefficient when we deal with a large number of openings and/or faults. To overcome this problem, we propose an iterative method which allows us to significantly increase the number of boundary elements and consequently the number of panels and faults by making the calculation separately for each fault and for each opening. During the iteration procedure, stress corrections are done in fault planes, with respect to two criteria, in order to account for permanent traction or shear failure (Mohr-Coulomb criteria for shear failure).

The code FAULT3D is presented herein, together with the scheme of calculation. A model verification was carried out by comparison with a two dimensional model constructed with UDEC 2.0 and CESAR -LCPC 3.1. An application to the Provence coal mine location, is also presented to demonstrate the practical application of the new code.

Displacement discontinuities method (DDM)

The displacement discontinuity method is a variant of the boundary elements method for solving problems in solid mechanics in general [START_REF] Crouch | Computer simulation of mining in faulted ground[END_REF]. It is particularly appropriate for problems involving faults or joints, mining in tabular orebodies (which extend at most a few meters in one direction and hundreds to thousands in the two others). We only need to suppose, for both faults and panels in tabular orebodies, that the boundary consists of two parts very close proximity to each other. The displacement discontinuity represents the relative displacement between the two sides of the boundary and has a physical meaning such as the closure of a mined excavation or the slip on a fault plan.

Let S represent a discontinuity surface within an elastic, homogenous, isotropic and infinite solid (Figure 1). X Q : field point From the general theory [START_REF] Linkov | Computer aided in analysis of stressed state and rockburst hazard in veins and coal seams influenced by faults[END_REF], it follows that the induced stresses at any field point Q are given by the integral : i, j, k = X,Y,Z [START_REF] Heib | Numerical modelling applied to dynamic phenomena -2 nd North American Rock Mechanics Symposium -June 1996 -Numerical modelling of ground burst problems -Workshop notes[END_REF] <j > "d(Q) : induced stresses tensor at the field point Q.

• O ® P sou» e r oint

AU k -

: k th k -Ul -Ul : kt h displacement discontinuity component at the integration point P. U + (IT) is a limit value of displacements when approaching S from the side for which the normal to the surface at the point P is an outward (inward) normal. R : influence coefficient of the integration point P on the field point Q. E : Young's modulus. v : Poisson's ratio.

If Q is tending to a point Q o of the surface S, we obtain in the limit of the surface :

where T ind (Qp)= On X\ Xi [START_REF] Crouch | Computer simulation of mining in faulted ground[END_REF] =0-i f (Qp).nj(Qo), a n is the T\ and Ti are the shear traction boundary at the point Qo, nj(Qo) is the j th component of the normal to the surface at the point Q o .

Let us write S as a sum of N elements I 5-I Each of them is called a displacement discontinuity element. The induced traction at the Q th element (Q = 1,.., N) is obtained from [START_REF] Crouch | Computer simulation of mining in faulted ground[END_REF] by : i, j, k = X,Y,Z AU P : is the constant displacement vector at the P th element.

If we know a method to evaluate the integrals over the elements, and if the displacement discontinuity vector at each element (3 components) is known we can find using relation [START_REF] Fotoohi | Nonlinear boundary element analysis of rock mass with discontinuities[END_REF] the traction boundary induced vector at any point in the surface.

In the next section, we explain how to use this method to analyse mining rock mechanics problems involving major discontinuities or faults

Use of DDM in mining rock mechanics problems

The key concept underlying the DDM is considering the tabular deposit or the fault plane as a discontinuity in a homogeneous, isotropic, linear elastic rock mass. DDM, unlike other numerical methods (finite elements for example) only requires the subdivision of the problem boundary (tabular deposit or fault plane) into elements instead of meshing the whole rock mass. Thus, it requires much less computational effort than any other equivalent differential method.

The boundary elements are considered in two different ways depending if they belong to a mined area or a fault plan. In the first case, we assume that one face corresponds to the roof and the second to the floor of the excavation. In the second case, we suppose that the opposite face of a discontinuity are connected with springs whose normal and shear stiffnesses are known.

Fault plane

Tabular orebody

Figure 2 : Typical mining problem solved by the DDM traction the normal to the surface.

Mined area

We can say that the total stress in an element, let us say q, is obtained as the sum of the initial stresses and the induced stresses by all the elements in the mined area. For these elements, the total stresses are zero. We have,

T T <"=T"+T h "'=(i -mo) ; = X,Y,Z [4]
This leads us to a system of 3*N equations with 3*N unknowns (three components of the displacement discontinuity vector for each element). If we suppose that the initial traction vector (3] components) in each element is known and have a method to evaluate the integrals over the elements, we can use an appropriate method to solve the system [START_REF] Sylla | Prévision du comportement des terrains sous l'influence des failles par des méthodes numériques et microsismiques -The 28[END_REF] to find the displacement discontinuity at each element. Using this solution, we can find via the integral [START_REF] Heib | Numerical modelling applied to dynamic phenomena -2 nd North American Rock Mechanics Symposium -June 1996 -Numerical modelling of ground burst problems -Workshop notes[END_REF] the complete induced stress tensor at any field point around the mine opening.

Fault plane

In this case, it was assumed that the fault does not slip until the excavation is made in its vicinity. For modelling purpose, a fault can be considered to be filled by a compressible material. Thus, the opposite faces of each fault elemental displacement discontinuity can be modelled to be connected by a spring chosen to be representative of the properties of the filling material (Goodman, 1976). The value of the displacement discontinuity components at each element, therefore, will be related to the normal and the shear stresses acting on it. The normal and the shear induced boundary tractions ran be related by the stiffness matrix to displacement discontinuity by (K. Fotoohi, 1993) [START_REF] Fotoohi | Nonlinear boundary element analysis of rock mass with discontinuities[END_REF]:

find = find KnO0 1 AUx o o^J[ A f/ z [5]
From (3) we can write, K,, 0 0 11 f

QK.O lAUr 0 0 K,l^U z 0=. \Rj t (Q,P)nj{Q)dS * SPO.SU Where, 0 y = 1 if i=j ^ lerwise. S k = K n if k=X and S k = K s if k = Y or Z. Q=1 N.
It simply corresponds to add K n and K s to the diagonal terms of the fault's influence matrix. The left side of the system is zero in respect with the assumption that there is no fault plan displacement before the excavation.

Consideration of the non-linear behaviour of the faults should lead us to more realistic results. For this purpose, two criteria have been chosen. These are presented below.

Tensile strength failure : cr n < R,

The normal stress o n on DD fault element can not exceed its tensile strength (R t ). When this condition is reached in DD element, the fault is considered open and acts exactly like an element on the boundary of an opening (Total stresses are zero).

Shear strength failure : T < C + a n tg (<p )

The Mohr-Coulomb condition is chosen to model the shear behaviour of the fault. The meaning of this constraint is that the total shear stress cannot exceed a certain amount given by the relation C+ontg((p), where C is the cohesion of the fault, <p its friction angle and o n the normal stress acting on the fault surface. If an element fails according to this criterion, the fault is considered in permanent slip and the excess shear stress is taken as the new boundary condition in the relation [7] instead of zero.

Verification of the code

In this section a verification of FAULT3D code is presented by comparison with the distinct element code UDEC 2.0 and the finite element code CESAR-LCPC 3.1. Two models involving the mining of a single panel, are presented here (figure 4). The width of the panel is fixed to 200 m and is 1000 m long, thus to justify plane stress analysis. The seam depth is 1000 m. In the first case, we consider two faults parallel to the panel, along its length, and in the second one, only one of the two faults is parallel to the panel. The results of a two dimensional modelling depend, in the latter, on where the cross section is chosen.

Rock mass properties :

The elastic rock parameters are : Young's modulus 10,000 MPa, Poisson's ratio 0.1 and unit weight 0.025 MN/m 3 .

Faults parameters :

The above data present the geomecanic parameters of the two faults. In situ stresses

The initial stresses* correspond to the overall virgin stress field. It is assumed that they are not modified by the tectonic activity nor by faulting. At this depth (1000 m) and with a unit weight equal to 0.025 MN/m 3 , the in situ stresses are : a v = CFE-W = ON-S = 25 MPa

Analysis of the results

Figure 5 shows that only the bottom part of the right fault is on shear failure. This part lies in the unloaded region of the panel. This failure corresponds to a decrease of the normal stress (combined to a small friction angle on the right fault) causing a lower failure limit in the Mohr-Coulomb criteria. Looking at figures 7 and 8, we can also observe that all the fault elements at the bottom of the panel are also in shear failure. All these three codes did not show opening elements in the faults.

If we consider the figure 6, which corresponds to the case 2, in a two dimensional modelling the results will be different from one cross section place to another.

Here we have supplementary information in the axial direction.

FAULT L FAULT R

Shear failure zone We can conclude, from this case, that the results given by FAULT3D are correct. Even if the calculations are in 3D, they require, here, a very low computational effort in terms of time of calculation and input data for a three-dimensional model.

Case study : Application to a panel of Provence colliery in a faulted zone

The Provence colliery is situated in southern of France between the cities of Marseille and Aix-en-Provence.

The archives found indicates that the mining began as early as the Middle Ages. Only one of the eight coal seams is presently extracted at a the depth of 1250 m. The dip of this seam is small (about 10"). The roof is quite competent which assures its relative stability. This rigidity is however the source of rockbursts caused by sudden release of stored energy due to the previous stresses and the induced stresses by the excavation [START_REF] Heib | Numerical modelling applied to dynamic phenomena -2 nd North American Rock Mechanics Symposium -June 1996 -Numerical modelling of ground burst problems -Workshop notes[END_REF]. The other important characteristic is the presence of many faults around the mined area. For many years, numerical modelling, combined with seismic events, have been used to help predict and alleviate rockbursts or collapses due to fault slip burst [START_REF] Sylla | Prévision du comportement des terrains sous l'influence des failles par des méthodes numériques et microsismiques -The 28[END_REF] 4

.1 final panel length forecast in a faulted rock mass

This study concerns the extraction of a new panel, railed T05 (figure 9). The final decision to stop the mining of this panel will depend on the behaviour of the mine structure when it comes in the vicinity of a faulted area. It was decided to make this study using FAULT3D to determine the final length of the panel which would ensure safety and profitability. The major problem is to determined the fault's characteristics from their tracks observed in the seam. Geomecanical parameters are essentially determine by the filling material properties and the fault's geometry will be estimated by orientations in the seam. This study also includes the old panels mined out before, which are 702 and T03.

Description of the model

The study concerns three lengths of the T05 panel (450 m, 500m and 550 m). It takes into account four faults observed during the excavation of the galleries V29 and V31 (Figure 9). Figure 2 shows the geometry of the model. The geomecanical parameters estimated for these faults are listed in the table 2. This configuration induces a shear failure zone, only on the fault C, corresponding to 1% of the total surface of the fault plane. This zone is situated in the roof of the seam, behind the panel face and is caused by a decrease of the stress acting to the normal of this fault plane (figure 11).

Case 2 : length of 500 m Figure ( 12) shows an increase of the failure zone from 1 % to 5% of the total surface of the fault C. This zone is still situated behind the face and about 100 m of the panel is concerned.

The failure zone extends above the seam until the limit of the fault. With the impossibility to know the dimensions of the fault, the only way to determine the exact failure zone is to increase the extent of the fault to incorporate the entire failure zone. That was done by increasing the extent of the fault by 50 m. This led to a new failure zone corresponding to 8 % of the total surface.

Case 3 : length of 550 m

The progress of the panel toward the fault C induces the increase of the failure zone size to a new value corresponding to 11% of the fault surface (figure 13).

In the 3 cases there is no failure in the faults A, B and D.

It appears from this study that a length of the panel less than 500 m should not affect the fault stability. For an extraction length greater than 500 m, the propagation of fault C failure affects progressively the seam when the mining face approach the fault. Thus, this situation could trigger seismic activities leading to rockburst.

back analysis of observed phenomena during the T05 panel excavation

Fault slip failure depends strongly on the geometry and geomecanical parameters of the fault. We mentioned that all these parameters were very difficult to estimate. During the mining of this panel, four seismic events were recorded in the side wall of V31 gallery. Theses events were caused either by fault D or by the failure of the rock mass. However, the stresses induced by the progress of the panel are not great enough to cause these events. In the next section we will study the influence of fault D by modifying, essentially, its dip and its azimuth. We can see that the maximum of Sf is obtained when P is low and becomes minimum when p is near 90°.

We represent in figures 15 and 16, the modifications observed when 3 is equal to 75° and when its equal to 15°. In the first case less than 1% of the surface is in shear failure and in the second case 8% is in shear failure and 1 % in traction failure. The shear failure is caused by a decreasing of the normal stress to the fault and the traction failure by the stress concentration at the limit of the panel.

Figure 17 and 18 show a decrease of the failure zone when a is modified from 6° to 10°. Values of a greater than 10° reduce the size of Sf and consequently the number of elements prone to failure.

As a conclusion, we can say that in some cases the fault D can be at the origin of recorded events in the gallery V31. Knowing the mechanism of these events we can calibrate the fault parameters in order to retrieve the same failure mechanism on the fault and make our model more realistic for the forecasts modelling. While the geometry of the panels is known, the major challenge is to characterise the faults and their distribution and their orientations. Also, it is important to estimate the parameters that govern their plane displacements. We show that using other methods, seismic events recorded for example, to characterise the media behaviour can permit us, by back analysis, to adjust our model (faults parameters and geometry) for future modelling.

Figure 1 :

 1 Figure 1 : Definition of the displacement discontinuity

Figure 3 :

 3 Figure 3 : flowchart of calculation in FAULT3D

Figure 4 :

 4 Figure 4 : Definition of model verification code

Figure 5 :

 5 Figure 5 : Case 1 by FAULT3D Figure 6 : Case 2 by FAULT3D

Figure 7 :-

 7 Figure 7 : Case 1 by UDEC

FigureCase 1 :

 1 Figure 9 : map of T05 panel

Figure 11 :Figure 14 :

 1114 Figure 11 :panel length = 450 m Failure area

FigureA 3 -

 3 Figure 15 : a = 6, 0 = 75

Acknowledgement

It was helpful for us to get exploitation information from the Provence colliery engineers. We would like to thank them for their comments and useful discussions. This work has been financially supported by an industrial research from the Provence colliery; the authors are grateful for their support.


O Deck

M Ai

F Homand

EFFECT OF SUBSIDENCE ON STRUCTURES BUILT OVER OLD MINING WORKS IN FRENCH IRON BASIN

Ail mining works can generate strata movements and lead sometimes to serious consequences to structures. This paper presents an original analysis of mining subsidence effects on surface structures. We present in the first part of this paper the results of a numerical model, which allows a global understanding of the structure cracking process. A real subsidence profile was applied at the bottom boundary of a finite element model and soil-structure interaction was studied by introducing a simplified structure at the top of the numerical model. A critical analysis was then carried out with the previous parameters which were supposed to be responsible for rupture. The second part is about the built of a better numerical model. This goal led us to make a synthesis of different countries formulas used to estimate the subsidence parameters such as horizontal strains.

RESUME :

Toute exploitation minière peut engendrer des mouvements du terrain et entraîner des conséquences parfois graves sur les structures. Cet article présente une étude originale des conséquences d'affaissements miniers sur des structures de surface.

Nous présentons dans un premier temps les résultats d'une étude numérique qui a permis une compréhension globale du processus de dégradation. Un profil d'affaissement réel a été appliqué à la base du modèle et l'interaction sol-structure a été étudiée en introduisant une structure simplifiée à la surface du modèle. Nous procédé à une interprétation critique des paramètres qui sont supposés être à l'origine des dégradations. La seconde partie traite de l'élaboration d'un modèle numérique amélioré. Cet objectif nous a conduit à réaliser une synthèse de formules provenant de différents pays, permettant l'estimation des paramètres d'affaissements comme les déformations horizontales.

INTRODUCTION

Many cases of mining subsidence exist around the world. In Lorraine (France), some old iron mines subsided some years ago (Auboué (1996), Moutiers (1997) Roncourt (1999)) so that a lot of houses cracked. Consequences on the population were as much economical than psychological. These mines were exploited by the traditional rooms and pillars method. Subsidence occurs because of the failure of old pillars in areas in which pillars were not removed. There were more than a hundred meters deep and the recovery ratio was ranged between 0,36 and 0,55 in the three latest subsidences. Unfortunately many other villages stand upon old underground works so that the question of consequences of mining subsidence on structure and prediction of damages is asked. An empirical approach is usually used to predict damages according to horizontal strain and the length of the structure (NCB, 1975) i Tensile zone ;; Compressive zone j Horizontal strain : e \ azq:: In order to understand the phenomena of degradation, we decided to use a numerical modelling program. It allowed us to observe (he soil-structure interaction during the mining subsidence process. These numerical results reinforced us to improve the model in order to improve the numerical observations efficiency. It is used to consider the horizontal strain on the surface (Figure 1, Aissaoui 1999) to be responsible for damages [START_REF] Kratzsch | Mining Subsidence Engineering[END_REF]. The strain intensity is, hence, very important and seems to be an interesting parameter to compare results of the modelling with field observations. This explains our decision to list formulas of maximum strain estimation from other countries.

With these formulas and an analysis about the effectiveness of the numerical model r we propose aUhe-end of this paper a simply realistic model to solicit different kind of structures.

NUMERICAL STUDY

PRESENTATION OF THE MODEL

Our topic is to find a simple model ,that allows to observe consequences of soil strains, displacements and curvature on simple surface structure. Because of soil heterogeneities, it is impossible to predict subsidence profile which can occur in a real mining situation with a finite element software. So the purpose of this part is the presentation of an elastic, two dimensional isotropic model which allows to observe the soil structure interaction. The finite elements software "CESAR-LCPC 3.2.4" is used. The model consists of a long soil strip (500 meters long, 50 meters high) with a structure on its top (Figure 2). The structure is modelled by a short homogeneous strip (1 meter high and 25 meters long).

Numerical problem is solved with plane strains assumption.

So that structure is a 25 m large raft with unlimited proof. About the model r we agree-it-to-be really simplified and not completely realistic. For example, the soil stiffness is probably not the same along the depth. The foundation can not have the same behaviour as a real building (with foundations, beams, posts, windows) and its regular loading is not realistic because in a house, loading is concentrated under posts. So results must be compared with each other and numerical values do not have to be used alone.

RESULTS PRESENTATION

We first present results for the case of a twenty five meters long foundation with a Young's modulus Ef=20000 MPa, loaded by QF=0.15 MPa, located in the compression area (Figure 3). The first column presents results for a foundation without subsidence and the second one shows the interaction when subsidence occurs. Curves show parameter variations in the soil along a section under the structure (0.1 meter deep).

About vertical stresses, we observe the consequences of curvature. Compressive area is characterised by a concave geometry which leads to biggest compressive values at the edges of the raft and smallest at the medium. Horizontal stresses are very different between the initial situation of the raft and after subsidence. We can see no change just under the structure, but a really increase at the edges. This may be explained by the soil-structure interaction. Raft confines soil so that all compressive strains that can not occur under it, occur at the edges. We can see a very important increase of shear stresses, more than 600 %, which is responsible of foundation slipping at the edges. We can argue with such results about the importance of curvature in relation to horizontal stresses. Curvature solicitation seems to be secondary as the vertical stresses variation associated is very low (less than 15 %) in comparison with horizontal stresses and shear stresses variation (more than 500 %).

Such results are corroborated by literature assessments. [START_REF] Neuhaus | de la construction des maisons d'habitation en zones d'affaissements miniers[END_REF], Peng and al. (1981), Yokel and ai (1982), Bell (1988) After this presentation, we can discuss all results for the parametric study. Two statistical analysis have been made to evaluate relative importance of each parameter. The first one consist in a variance analysis, and the second one in an experimental design. We just set out results about relative importance of foundation characteristics to each other which is estimated by computing the stress variation rate when one parameter changes. The stress variation rate is more efficient than absolute change, because it may be compared to the safe coefficients used for the building construction.

Foundation without subsidence

We remind all the values for the parametric study in table 1. Shear stresses in the soil and solicitations in the foundation increase with the length. So length is very harmful. In compressive area, horizontal stresses at the edges of the structure will have very disastrous consequences. The Young's modulus seems not to have influence upon stresses. The reason is probably the big difference between the stiffness of soil and structure. In all cases, structure is dreadfully stiff in relation to the soil. But we notice that stresses in the foundation increase with the structure stiffness. The structure loading appears to be benefit because higher the loading is, smaller the variation rate is.

In conclusion, this study allows a better understanding of the phenomena. We have now to improve the geometry of the

HORIZONTAL STRAIN INVENTORY

The maximum horizontal strain evaluation is very important because it is supposed to be one of the more important parameters to explain damages (NCB 1975). So the comparison between numerical values and empirical values is a good way to check numerical strains first estimated.

Figure 4 shows values of proportional coefficient "K" to estimate maximum horizontal strain. This coefficient "K" must be multiply by a mining ratio "o/H" or "Vymax/H" in which "o" is the opening of the mining works, "Vymax" is the maximum vertical surface subsidence and "H" the mining depth (Figure 1). 1

-i 1-H 1 1 1 1 i 1 1 1 1 1 1 1 1-> K 1 1 1 (tension) 1 (compression) (tension) # lz (tension) 2 (tension) ' # 2 (compression)
e=K* o/H (gray) e=K * Vymax/H (black)

.

(compression)

. 2 (compression) .1

(tension) (compression) (tension) (compression)

# 5 (compression) . 5(tension)

Figure 4 : coefficient of proportionality for the maximum horizontal stress estimation.

If we compare these formulas with results of the first numerical model presented in the second section, we observe a really underestimated horizontal strain. Indeed, the first model was built with a real subsidence profile with a 150 meters deep mining work and a 1.2 m subsidence at the centre. Maximum strain should be about 1.2/150 = 8.10' a [eq. 1 with K=1]. Although we computed a maximum tensile stress of 0.15 MPa and a maximum compression stress of 0.27 MPa. Soil stiffness is of 100 MPa which leads to a maximum horizontal strain of about 1.5 10" 3 in tensile area and 2.7 10 3 in compressive area. The difference between both the results may be explained by the small model size which did not contained the mining works. In order to study soil-structure interaction we decided to use an improved numerical model which leads to a good strain rough estimate.

NUMERICAL MODEL IMPROVEMENT

We first want to estimate the errors due to the small height and the bottom boundary displacement condition of the previous numerical model. We decided to make a global numerical model which contained the underground mining works. From the results of this global model, a small model which the geometry was exactly the same than in the second chapter is solicited. By comparison of both results, it is possible to estimate the errors made when we try to reproduce real mining effects on the surface in the same way than in chapter two (that is to say with a small model, solicited by a subsidence profile at its bottom boundary). Maximum horizontal strain grows up to average 3.10 3 when "theoretical" value according to the horizontal strain formulas [eq. 1, with K=1], is 1/150 = 6.5.10 3 .

Results comparison e, surface

Vy, surface In fact, the two models react like a beam. If we impose a vertical displacement, horizontal strains at the top boundary will increase with the high of the model, indeed, a fix vertical displacement leads to a fix curvature. And we know that for a beam :

Mf.h/2 [2] a = J ' [3] Mf = E.I IT [4] a -E.e
With "a" normal stress at the top of the beam, "h" height, "I" moment of inertia, "E" Young's modulus, "1/R" curvature and "Mf" bending moment. Thus :

h {$] £ = 2.R
So when the height of the model increases, strains will grow up too. Then it is possible to find the optimal model size length "L" and height "h" with which numerical and real horizontal strains at the surface are the same. We saw that the height increase does not really change the vertical displacement profile, but the higher "h" is and the greater the computed strains will be. "L" must be find in order that the lateral boundaries do not disturb model results.

  w

Figure 1 :

 1 Figure 1 : Horizontal strain profile.

Figure 2 :

 2 Figure 2 : Numerical model presentation (Lf : foundation length, Qf : foundation loading, Et : foundation Young's modulus, E s : soil Young's modulus).

Figure 3 :

 3 Figure 3 : Stress variation along a 0.1 meter deep section in a compression area ; Ef = 20000 MPa, Lt = 25 m, Qi = 0,15 MPa.

  structure and to use an elasto-plastic criterion for the soil in order to avoid tensile stresses in the soil which are not realistic. The third improvement is about the bottom boundary condition from which the surface horizontal strains depend. In next section, we present our methodology to obtain a better model at the top of which it will be possible to study different kind of structures subjected to horizontal strains with vertical displacements.

  range value is refered by a number which allows to know the author. Authors references are the followings : (1) Wagner and Schuman 1985, (2) Kratzsch 1983, (3) Proust 1961, (4) Whittaker and Reddish 1989, (5) Yokel and al.1981. These coefficients result from coal mines in France, Great Britain, Poland, United States and Germany. We can see a big concentration around the unit value.

  Global model size is 200 m high and 570 m long. Mining works are 150 m deep, 200 m long and we simulate a 1 meter subsidence. To obtain such a vertical displacement at the centre of the subsidence profile, we impose a unit vertical displacement at the top of the underground excavation (Figure 5). This figure clearly explains our methodology. Vertical displacements computed at the top of the global model are used for the bottom boundary condition of the small model. Results are shown in Figure 6. We can see that we make strong errors when we only use a small model. Vertical displacements are correctly transmitted between the two models, but horizontal strains are divided by about 2.

Figure 5 :

 5 Figure 5 : Comparison between global and small models.

Figure 6 :

 6 Figure 6 : Vertical displacement and horizontal strain at the surface of the two models.

  from visual and practical observations allow to estimate structures damages with only few parameters, but these are not satisfactory. It is why we made numerical models which allow a better understanding of the structure cracking process. A first model allows to estimate the importance of structure length, stiffness and loading to each other. We showed that it was possible to use a simply realistic model to study the soil-structure interaction. In the future we think to improve again and valid this first study in order to use the model with different kinds of post, beam structures and to compare numerical results with those of a statistical study about effects of a real mining subsidence on houses. We hope to be able to improve empirical methods which are usually used to estimate damages. BIBLIOGRAPHY AISSAOUI K. (1999).-Amelioration de la prevision des affaissements dans les mines à l'aide des approches empiriques, numériques et analytiques. Thèse de doctorat INPL (Nancy). BELL F.G. (1988).-The History and Techniques of Coal Mining and the Associated Effects and Influence on Construction. In Bulletin of the Association of Engineering

Table 1 :

 1 Parametric values of the foundation characteristics.