Aurelie Martin
email: aurelie.martin@inrs.fr

Frederic Yves Bois

Francis Pierre

Pascal Wild

Aurélie Martin

Occupational Exposure to Cobalt: a Population Toxicokinetic Modeling Approach Validated by Field Results Challenges the Bei® for Urinary Cobalt

Keywords: biological monitoring, occupational health practice, mathematical models

The objective of this study is to model the urinary toxicokinetic of cobalt exposure based on of 507 urine samples from 16 workers, followed up for one week, and 108 related atmospheric cobalt measurements to determine an optimal urinary cobalt sampling strategy at work and a corresponding urinary exposure threshold (UET). These data have been used to calibrate a population toxicokinetic model, taking into account both the measurement uncertainty and intra-and inter-individual variability. Using the calibrated model, urinary sampling sensitivity and specificity performance in detecting exposure above the 20 μg/m 3 threshold limit value (TLV-TWA) has been applied to identifying an optimal urine sampling time. The UET value is obtained by minimizing misclassification rates in workplace exposures below or above the TLV . Total atmospheric cobalt concentrations are in the 5 -144 µg/m 3 range, and total urinary cobalt concentrations are in the 0.5 -88 µg/g creatinine. A two-compartment toxicokinetic model best described urinary elimination. Terminal elimination half-time from the central compartment is 10.0 hours (95% confidence interval (8.3 -12.3)). The optimal urinary sampling time has been identified as 3 hours before the end of shift at the end of workweek. If we assume that misclassification errors are of equal cost, the UET associated with the TLV of 20 µg/m 3 is 5 μg/L which is lower than the ACGIH recommended BEI® of 15 μg/L.

Introduction

Cobalt (Co) enters into the composition of hard metal alloys used in electrical, aeronautical or car industries as a binder for metallic carbides (tungsten, silicon, vanadium). Cobalt enhances the resistance of these alloys to high temperatures and corrosion. Cobalt is also used in steel production, glass and ceramic industries, and many other industrial applications. All personnel working in these industries are potentially exposed to cobalt or its compounds, both of which can induce chronic intoxication leading to respiratory, thyroid, cutaneous, cardiac and carcinogenic effects [START_REF] Lison | Experimental research into the pathogenesis of cobalt/hard metal lung disease[END_REF][START_REF] Demedts | Respiratory diseases from hard metal or cobalt exposure .Solving the enigma[END_REF][START_REF] Shirakawa | Occupational asthma from cobalt sensitivity in cobalt workers exposed to hard metal dust[END_REF][START_REF] Lauwerys | Industrial toxicology and occupational poisoning[END_REF][START_REF] Moulin | Lung cancer risk in hard metal workers[END_REF] under some circumstances.

The American Conference of Governmental Industrial Hygienists (ACGIH) has defined two exposure indices applicable to cobalt occupational exposures (6): an atmospheric threshold limit value time weighted average (TLV-TWA) of 20 µg/m 3 and a biological exposure index (BEI®) of 15µg/L total cobalt concentration in urine collected at the end of shift at the end of the last day of the workweek. The ACGIH has justified this BEI® through four studies [START_REF] Scansetti | Urinary cobalt as a measure of exposure in the hard metal industry[END_REF][START_REF] Lison | Biological monitoring of workers exposed to cobalt metal, salt, oxides, and hard metal dust[END_REF][START_REF] Ichikawa | Biological monitoring of cobalt , based on cobalt concentrations in blood and urine[END_REF][START_REF] Swennen | Epidemiological survey of workers exposed to cobalt oxides, cobalt salts, and cobalt metal[END_REF] giving linear regression equations between cobalt concentration in ambient air and in urine obtained at the end of shift at the end of the workweek. In each of those studies, the regression equation was used to estimate cobalt levels in urine at 20 µg/m 3 exposure, assuming the TLV value to be reliable. No data on exposure or urinary cobalt were available before the last shift, except in one study [START_REF] Scansetti | Urinary cobalt as a measure of exposure in the hard metal industry[END_REF]. Thus the potential effect of cobalt exposure during the preceding days could not be assessed. Furthermore, in the absence of accurate information on cobalt excretion kinetics and of inter-individual variability, it is unclear whether urinary collection at the end of shift at the end of workweek is indeed the best sampling time. The BEI® value obtained from a linear regression equation is itself certainly not optimal because it ignores both the time-dependency of relationship between airborne exposure and urinary cobalt and the very large inter-subject variability. This paper presents a toxicokinetic approach for determining an optimized urinary sampling strategy for cobalt and a consistent way of deriving an urinary exposure threshold (UET), corresponding to the ACGIH BEI®, for use in occupational hygiene. A population toxicokinetic model of cobalt intake and urinary excretion was developed from data on a series of cobalt-exposed subjects followed up during a period of one week. Atmospheric and urinary cobalt measurements were simultaneously collected for each subject. The model was then applied predictively to identify an optimal sampling strategy and back-calculate a UET corresponding to the current value of the TLV, which is assumed to provide valid protection for exposed workers.

Population and Methods

Population

A total of 16 male subjects exposed to cobalt dust were recruited in two plants producing tungsten carbide cutting tools. They were followed-up during one workweek respectively in 1988 and 1993.

Atmospheric and urinary sampling strategy

Airborne exposure to cobalt was measured by sampling the inhalable dust fraction in each individual breathing zone for all work shifts during the study week. The workweek comprised four, five or six days, depending on the plant and the subjects. For most subjects, exposure was measured for each half-shift (mean duration of half-shift = 3.5 hours), and for the others, exposure was measured for the total shift (mean of total shift = 7.5 hours).

Airborne sampling was performed in compliance with French AFNOR standard NF X 43457. Aerosol samples were collected on WHATMAN® QM-A quartz fibre filters fitted in closed configuration inside a MILLIPORE 37 mm cassette designed to collect the aerosol inhalable fraction by means of GILIAN® HFS113 (1 L/min flow) portable sampling pumps.

The personal pump flow rate was carefully checked at the beginning and end of each sampling sequence. The flow rate was also regularly checked at intermediate times.

All urine voiding of the study week were collected in parallel for each subject. Each subject was informed of the study aims and conditions, and agreed to collect urine based on a controlled sampling procedure: Subjects were requested to collect urine samples whilst avoiding contamination, without working clothes, with clean hands and, if possible, after a shower at the end of shift. The urine was processed, packaged and frozen on site. All equipment cups and tubes in contact with the urine for analysis were previously washed with a hot detergent, rinsed and immersed for 48 hours in a 10% nitric acid solution , rinsed with ultra-pure water and dried in a oven at 50°C. All these operations as well as storage were performed in airtight polyethylene containers, in which all devices were protected from dust and contact with operator fingers. Each urine sample was analyzed separately. For 12 subjects, all voiding of the following weekend were also obtained and analyzed.

Analytical methods

Air sample analysis

Total cobalt content of the filters was analyzed by two laboratories with a method developed by one of them [START_REF] Demange | Aerosol evaluation difficulties due to particle deposition on filter holder inner walls[END_REF]. These laboratories belong to the French inter-laboratory ring trial network for occupational hygiene. Particles deposited on each filter, including dust on the filter holder inner walls, were dissolved in mixture of hydrofluoric (2 mL, 40%) and nitric (3 mL, 68%) acids. Cobalt was measured at a wavelength of 240.7 nm by atomic absorption spectrophotometry, using a flame technique based on a 10 µg/L detection limit.

Urine sample analysis

Total urinary cobalt concentration was evaluated by electrothermal atomic absorption spectrometry using a Perkin Elmer 3030/Zeeman instrument based on a method [START_REF] Baruthio | Cobalt determination in serum and urine by electrothermal atomic absorption spectrometry[END_REF] involving a detection limit of 2 nmol/L ≈ 0.1µg/L in INRS internal quality system. Urinary creatinine was determined by colorimetry (Roche 3667 kit) using a Cobas Mira S Plus (Roche Diagnostic System).

Statistical methods

A number of urine sampling strategies were compared on the basis of their efficiency at detecting an atmospheric cobalt concentration higher than a limit value. Optimal sampling was determined based on results of a multilevel (population) toxicokinetic data modeling [START_REF] Bois | Population toxicokinetics of tetrachloroethylene[END_REF], which estimates jointly the toxicokinetic parameters of each subject (subject level) and the inter-worker distribution of these parameters (population level).

Population toxicokinetic model

Figure 1 shows a conceptual graph of this model.

At the subject level, a toxicokinetic approach was used to model the urine-excreted quantity of cobalt in each voiding for each subject, as a function of the time-dependent atmospheric cobalt concentration in the worker's breathing zone, and of the urine sampling time.

All urine voiding were collected and analyzed separately. The absolute quantity of cobalt (in µg) excreted at each voiding was modeled as the dependent variable. This quantity was obtained by multiplying the cobalt concentration by the corresponding volume of excreted urine. Thus, this process overcame the need to use cobalt expressed per g creatinine.

For each subject, the time course of cobalt quantity in urine (Q u ) was modeled as the output of a deterministic two-compartment model. This model included a central compartment and a peripheral one, both without particular physiological interpretation. We checked that a simpler 1-compartment model would not correctly predict the data (results not shown). We also checked different parameterizations using various distribution shapes (see the Discussion section). The toxicokinetic model for each subject i had four transfer parameters θ i = {K in , K r , K s , K e } i characterizing transfers between compartments (the population distributions for K r i and K s i are not displayed on Figure 1 for clarity). The mathematical model and the precise meaning of these parameters are provided in online appendix 1. The toxicokinetic model can be viewed as a function relating cobalt urinary excretion to air exposure, time, and toxicokinetic parameters for a given subject:

Q u = f(C in , t, θ i ) (1)
At the population level, subjects were assumed to differ randomly from each other.

The 4-component parameter vectors θ i characterizing each subject were assumed to be lognormally distributed in the population:

Log(θ i ) ~ Normal(μ, Σ) (2)
in which μ and Σ were respectively the population mean and variance (themselves vectors of four elements each, the variance measuring inter-individual variability) in a logarithmic scale.

At the measurement level, both cobalt concentration in air C in and urine cobalt quantity Q u were measured at finite accuracy. A measurement error model, assumed applicable to all subjects, was therefore set up to account for uncertainties affecting those data. We assumed that C in was measured with a multiplicative log-normal error with GSD 1.5

(implying a 95% chance that the measured value was between 0.5 and twice the true value) (see the discussion for a justification of that model). The analytical measurement errors around Q u were assumed to follow a normal distribution with mean zero and a standard deviation, ξ, modeled as the sum of a constant error term SD min and a term proportional to Q u (eq. 3):

CoV Q SD u ⋅ + = min ξ (3)
When Q u is close to zero, ξ is at least equal to SD min and, for high values of Q u , that equation approximates a relative error model. The proportionality term (CoV) can be interpreted as a coefficient of variation.

All the parameters of the above models (θ i , μ, Σ, CoV, SD min ) were estimated in a Bayesian framework. Because for most parameters little prior information was available (mostly that they are non-negative and their order of magnitude), their prior distributions were chosen so that they had little influence on the final result (see online appendix 2). These prior distributions were then updated on the basis of the data to obtain posterior distributions, which were the Bayesian equivalent of estimated parameter and confidence intervals. Updating required calibration of the entire statistical toxicokinetic model with urinary and atmospheric cobalt concentrations data measured for the 16 subjects. Markov Chain Monte-Carlo (MCMC) methods (14-15) were applied using MCSim software (Free Software Foundation, Boston, MA, "http://www.gnu.org/software/mcsim/"). Two parallel Monte-Carlo Markov Chains with different starting points were run. After 10000 iterations, the two chains mixed well and converged, according to the R criterion of Gelman and Rubin (1992) [START_REF] Gelman | Inference from iterative simulation using multiple sequences (with discussion)[END_REF]. The following 50000 iterations of the two chains were used for identifying the posterior distribution. Model fit was checked and its predictive properties were tested by crossvalidation (details on priors and model checks given in online appendix 2).

Determining an optimal sampling strategy

To determine the optimal strategy, we randomly sampled the toxicokinetic parameters for 200 virtual individuals from the posterior population distributions obtained as described above. Each of these virtual individuals was exposed to a cobalt given dose for one workweek (four 8-hour shifts). This dose was assumed to be equal for all shifts. For each subject, a series of 28 atmospheric cobalt exposures was set, with exposure concentrations between 5 μg/m 3 and 30 μg/m 3 For any USS, the excretion rate of cobalt in urine was simulated every hour and the cobalt content of voiding, assumed to take place every three hours (close to the observed mean duration), were computed. A random error was added to that quantity based on the above urinary measurement error model (eq. 3). Predicted cobalt quantities were transformed into urinary cobalt concentrations by dividing them by 237 mL, the mean observed volume among the study subjects. That procedure takes into account the dependence of cobalt concentrations (inside the body and in urine) on the toxicokinetic parameters of each individual, and on the time-varying exposure (cobalt concentrations increase nonlinearly during exposure and decrease in the absence of exposure).

USS were compared according to their sensitivity and specificity in detecting an atmospheric value above a TLV of 20 μg/m 3 . For each USS, a receiver operating characteristic (ROC) curve, giving the sensitivity for detecting above TVL excursions as a function of (1 minus specificity), was constructed by varying the urinary exposure threshold (candidate UET) by steps of 0.25 μg/L between 0 and 9 μg/L and by steps of 1 μg/L between 10 and 20μg/L. For any given value of that threshold (say 8 μg/L), the sensitivity was computed as the number of simulated urinary concentration values exceeding 8 μg/L, divided by the number of simulations for atmospheric exposures equal or superior to the 20 μg/m 3 TLV (n=13x200=2600). Similarly, the specificity was computed as the number of simulated urinary concentration values below 8 μg/L, divided by the number of simulations for atmospheric exposures below the 20 μg/m 3 TLV (n=15x200=3000). The best USS was taken as the one with the largest area under the ROC curve [START_REF] Hanley | The meaning and use of the area under a receiver operating characteristic (ROC) curve[END_REF].

After determining the optimal USS, the final step was to establish an optimal UET value for comparing future urinary results. For each candidate UET, each of the 5600 values of the USS-specific urinary measurements was either well-classified (both the atmospheric or urinary are either below or above their respective limits), or falsely positive (i.e., with an urinary value above the UET with an atmospheric exposure below the TLV) or falsely negative (i.e., with an urinary value below the UET with an atmospheric exposure above the TLV). We denote by C FP the cost associated with a false positive classification, C FN the cost associated with a false negative classification and FP and FN, respectively, the number of false positives and false negatives in our sample of 5600 values. Denoting by C= C FP /(C FN + C FP ), the relative cost of false positives over all false classifications, each value of the UET was associated with a total cost, given by:

Total Cost= CxFP+(1-C)xFN

A large UET resulted in wrongly classifying high urinary exposures as acceptable (large FP and low FN) while a small UET resulted in wrongly classifying low exposures as unacceptable (large FN and low FP). Increasing C, i-e increasing C FP whilst keeping C FN constant, would lead to put more emphasis on the costs associated with wrongly deciding that the atmospheric exposure exceeds the TLV based on the urine sample, while decreasing C would put more emphasis on the protection of the worker. The optimal UET among the candidate UETs we considered (between 0 μg/L and 20 μg/L) , corresponded to a minimum total cost.

Results

Population characteristics

The atmospheric and urinary measurements for each of the 16 subjects studied are summarized in Table 1. The geometric mean of the atmospheric cobalt exposure measurements of plant A (subjects A to G) exceeded the TLV of 20 µg/m3 defined by the ACGIH in 3 out of 7 subjects. In plant B (subjects 1 to 9), the geometric means of the atmospheric cobalt measurements were even higher. Atmospheric cobalt exposures for each subject were usually followed by a peak in urinary cobalt within the next hours, which declined until the next exposure. Urinary cobalt excretion declined gradually over the week end.

Modeling results

The 2-compartment model gave a relatively close data fit, with model predictions following the described previously concentration time-course pattern. We chose the best model among the different ones we tested using different priors and distribution shapes.

Those made little difference, but as soon as we chose a two-compartment model, the model fit improved markedly and did not depend much on the distributional shapes and priors we used.

The overall correlation between observed data and model predictions is given in online appendix 2. Figure 2 illustrates in detail the fit for subject A. To provide a more familiar representation of the kinetics illustrated in Figure 2, all urinary measurements and predictions (urinary quantities in µg) were divided by time elapsed since the previous urine collection. This led to an expression of the urinary cobalt excretion flow in µg/h. Table 2 presents a summary of the posterior distributions of population and individual model parameters.

The population geometric mean of K in characterizing the quantity of inhaled cobalt entering the organism per unit time, was estimated at 60.4 L/h with a between-subject GSD of 2.0. Elimination parameter K e corresponds to an elimination half-time T e1/2 of approximately 10 hours, varying between 8 and 12 hours relatively constant for the studied subjects (GSD=1. [START_REF] Gelman | Inference from iterative simulation using multiple sequences (with discussion)[END_REF]). Half-times T r1/2 and T s1/2 corresponding to exchanges with the peripheral compartment, were estimated at 9 and 20 hours respectively, with a rather high interindividual variability. The coefficient of variation (CoV) for the urinary cobalt measurements was estimated at about 36% (CI 95%: 32% to 40%), with a baseline SD min of 0.03 µg/voiding (CI 95%: 0.002 to 0.08)

Optimal sampling strategy

Detailed comparisons of the various strategies, on the basis of ROC curves, are presented in the online appendix 3. The best sampling strategies were to collect urine sampled during the last 3 hours of the last shift of the week (USS D) or during the 3 hours following the end of the last shift of the week (USS E). For ease of sampling, urine sampling from the last 3 hours of the last weekly shift can be considered as the optimal sampling strategy. Figure 3 illustrates cost functions for two misclassification costs, C. The heavy line represents the cost associated with a false positive (FP, wrongly deciding that the exposure exceeds the TLV) and a false negative (FN, wrongly deciding that the atmospheric exposure is below the TLV). These were set to the same value, meaning that the user gave as much importance to an FN as to an FP (C = 0.5 and C FN =C FP ). The optimal UET value is 5 µg/L in that case. For the exposure values considered and given the chosen optimal USS, this UET yields an FP percentage of 40% and of an FN percentage of 23%. For the light line, the cost of FN was considered to be twice as high as that of FP meaning that the user gave twice as much importance to an FN as to an FP (C = 0.33 and C FN = 2C FP ). The optimal UET was then approximately 3.5 µg/L, yielding 58% FP and 12% FN.

Discussion

This paper presents an optimized strategy for urinary sampling in the workplace based on a population toxicokinetic modeling of the urinary cobalt excretion kinetics in 16 workers over a period of one week. We find that the optimal strategy was to sample urine during the last three hours of the workweek. In implementing this strategy, we have calculated a UET value corresponding to the TLV of 20 µg/m 3 depending on costs assigned with wrong decisions. A UET value of 5 µg/L is obtained when the cost of wrongly deciding that exposure is acceptable is considered equal to the cost of deciding that the exposure is unacceptable. Our results agree with international recommendations (6) with respect to the urinary sampling time (end of shift at the end of workweek). They also suggest that sampling after the end of shift at the end of workweek would not provide any additional benefit. However, the optimal biological value (UET) we obtained (5μg/L or 3.5 μg/L, depending on the error relative cost) is lower than the value (15 μg/L) recommended by the ACGIH. When applying the 15 μg/L value, we estimated 80% of false negatives (that is, workers with urinary cobalt concentrations below the UET during exposure atmospheric concentrations exceeding the TLV), and 3.5% of false positives (data not shown). The high percentage of false negatives shows that the ACGIH BEI® may not sufficiently protect the workers. The difference depends probably to a large extent on how correspondence between TLV and BEI® was established. The ACGIH have used only cross-sectional end-of-shift urinary data linearly regressed on the atmospheric exposure. Most notably, no individual kinetics were observed. Conversely, our results stem from modeling of actual individual data for workers followed-up over a full week. It should be noted (data not shown) that, if we apply the ACGIH strategy to our data with a simple linear regression, the corresponding value of the UET is 15 µg/L when we regress urinary concentration on the log atmospheric concentration in the last shift of the workweek, and is 12 µg/L when using atmospheric concentration on the natural scale. If we apply the ACGIH strategy to our data, we obtain an UET close to the BEI®. However, our data and method also allow to actually estimate misclassification rates and, allowing for them, yield a much lower UET. This is largely due to the fact that our approach takes inter-worker variability into account; the ACGIH approach does not do this explicitly. Thus our results represent a challenge not only to the BEI® for cobalt, but, more generally, to the ACGIH approach to deriving these BEI®s.

Our results may be discussed in relation to a number of issues. We did not consider varying volume of voiding when determining an optimal strategy and the corresponding UET. However, we do not believe that this is a serious limitation. What was measured, and used as data to calibrate the model, is the cobalt quantity in µg in each voiding. That quantity depends on the time between voiding (the voiding times were recorded and those actual times were used in input to the model). The urine volume also depends on time between voiding, but the quantity excreted between two voiding is conditionally independent of the urine volume given the voiding times. However, we used the mean volume observed among our subjects, when converting this quantity into a concentration in µg/L, for display purposes.. If we wanted to ensure protection for all workers, we would need to use a minimum urine volume and this would lead to even higher concentrations.

All urines of the week were collected and atmospheric exposure was measured for all the shifts worked. The range of cobalt urinary and atmospheric concentrations was very wide (geometric mean from 4.89 to 144.22 µg/m 3 for atmospheric data, and from 1.62 to 12.33 µg/g creatinine for urinary data). The range of applicability of our results should therefore be relatively wide too. Moreover, the total number of workers ( 16) studied is not very large, but for each worker all exposure and urinary follow-up is complete.

The model estimates at 36% the coefficient of variation of the urine analytical cobalt concentration measurement . This is higher than the laboratory value, supposedly to be 4% [START_REF] Baruthio | Cobalt determination in serum and urine by electrothermal atomic absorption spectrometry[END_REF]. We therefore confirm that the true accuracy of field studies is probably lower than pure laboratory uncertainty. However, the estimated 36% CV also includes other sources of errors (modeling error or intra-individual variability) and represents an upper limit in terms of analytical accuracy.

Atmospheric exposures measurements were collected over half shifts or complete shifts because instantaneous readings were not feasible. Therefore, in contrast to usual pharmacokinetic modeling, only mean exposures values were available in this case. This created a degree of uncertainty in the estimated parameters. Secondly, atmospheric measurements cannot be assumed to correspond exactly to the inhaled quantity. This was taken into account by assuming a log-normal measurement error with GSD 1.5 around the true inhaled concentration value. Other GSDs between 1.2 and 2 were tried, but the model's overall fit did not change much and that had no influence on the UET. Only one model is presented in this paper. However, alternative toxicokinetic models were considered and fitted, including one-compartment models and other parameterization of the measurement submodels. For instance, several different prior statistical distributions were tested. The model shown here was chosen because of its better fit to the data.

The parameters of our model are subject to only limited physiological interpretation.

The elimination half-time T e1/2 (equivalent to rate constant K e ) can be interpreted as the body elimination half-time of cobalt (central compartment). The parameters determining flow rates between the central compartment and other non-specified organs (lumped into a single peripheral compartment) cannot be easily interpreted. The estimated geometric mean for K in , the parameter determining the cobalt quantity inhaled into the organism per unit time, is 60.4 L/h. That is quite lower that the physiological respiratory flow rate, which is about 500 L/h for a man at rest. A possible explanation for this is that 88% (1-(60.4/500)=0.88) of the inhaled cobalt is exhaled or does not enter the body (i.e. is eliminated by the lungs into the gastro-intestinal tract and feces without absorption). It has been estimated that approximately 30% of cobalt inhaled as cobalt oxide can be absorbed [START_REF] Elinder | Cobalt. In Handbook on the Toxicology of Metals[END_REF].

It is interesting to note that, when a single compartment model was applied (data not shown), the estimated parameter K in was virtually identical to the one estimated in the twocompartment model. Elimination half-time from the one-compartment model was approximately 20 hours (data not shown), which is close to the values given by Lauwerys and Hoet (19) (p88) quoting Christensen et al [START_REF] Apostoli | Urinary cobalt excretion in short time occupational exposure to cobalt powders[END_REF] and Apostoli et al [START_REF] Christensen | A short-term cross-over study on oral administration of soluble and insoluble cobalt compounds: sex differences in biological levels[END_REF]. In the twocompartment model, the apparent elimination time depends on the rate-limiting exit from the peripheral compartment rather than on the elimination from the central compartment. But here also , the elimination half-life is estimated at approximately 20 hours (T s1/2 =20.38 hours) whilst the 2-compartment model fit to the data was far better than that of a one-compartment model.

The results given in this paper were based on a large number of data (507 urinary cobalt determinations and 108 atmospheric cobalt measurements), requiring very close cooperation of the study subjects and availability of research personnel for every working shifts during the week. Such a protocol is expensive and difficult to organize. Therefore, it should only be used if its results are expected to be interpretable and useful. An important constraint in this respect is the elimination half-time. If it is longer than a few days, this procedure involving atmospheric and urinary data for a whole week will not allow it to be accurately identified it with any precision, yet it would still condition the overall body burden. Thus for metals with a longer elimination half-time, this protocol would be unsuitable. A similar protocol with sparser sampling over a longer period of time would be more appropriate.

Optimal design methods could be applied to this issue (22). Lauwerys and Hoet [START_REF] Lauwerys | Cobalt. In Industrial chemical exposure[END_REF] suggest there is another component (possibly due to kidney or liver storage) of cobalt kinetic, which may persist for 2 years. The protocol applied in this study would naturally be incapable of identifying phase associated with such a compartment. However, for all practical purposes in occupational health, two compartments described the data and simulations closely enough (data not shown), showing that a dynamic steady-state is reached.

On the other hand, we should note a number of limitations on the possibility of extrapolating our study results. First, despite the wide range of exposure levels, all data were obtained in hard metal factories at which cobalt was always combined with tungsten. It cannot be assumed that the results would be exactly the same in other forms of cobalt exposure (coating, recycling, ceramics, polymers). However, the data do originate from two different plants, which contributes to our confidence in results representative of such this form of occupational exposure. It is known that the close-faced cassette sampler is slightly biased in relation to the ISO curve. However, the exposure is never overestimated so the low percentage of cobalt retained in the body cannot be explained by this bias. Furthermore, this study assumes that all the cobalt intake is via inhalation. It is probable that some intake is via ingestion, on which we have no information [START_REF] Lauwerys | Cobalt. In Industrial chemical exposure[END_REF]. However, the low percentage of cobalt retained in the body does not support the hypothesis of a major impact from a source of exposure other than the respiratory. A final limitation of our study is that no physiological data were available on the subjects (e.g., respiratory flows, smoking habits and so on).Therefore, individual characteristics could not be taken into account using, for example, a physiologically-based pharmacokinetic model. However, this may not be central issues, given the aim of our study (obtaining valid USS and corresponding UET). The corresponding coefficient of determination is 0.77 on the natural scale and 0.72 on the logarithmic scale. Although some of the measured values are below the corresponding predicted values, the overall fit seems adequate. Note that (with the possible exception of subject 9 which has one of the lowest atmospheric exposures) the various "outliers" correspond to different subjects.

E E E E E E E E E E E E E E E E E E E E E E E E E E E E E E E E F F F F F F F F F F F F F F F F F F F F F F F F F F F F F F F F F F F F G G G G G G G G G G G G G G G G G G G G G G G G G G G G G G G G G G G 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 2 2 2 2
For atmospheric measurements, one clearly sees the effect of rounding in the lower measurements. The overall fit is however adequate. The corresponding coefficient of determination is 0.82 on the natural scale and 0.88 on the logarithmic scale.

E E E E E E E F F F F F F F F G G G G G G G G 1 1 1 1 1 1 1 1 1 1 2 2 2

Comparison of ROC curves for various Urinary Sampling Strategy (USS)

As explained in the main text, the USS were compared on the basis of their sensitivity and specificity in detecting an atmospheric value above a TLV of 20 μg/m 3 . For each USS, a receiver operating characteristic (ROC) curve, giving the sensitivity for detecting above TVL excursions as a function of (1 minus specificity), was constructed by varying the urinary decision threshold by steps of 0.25 μg/L between 0 and 9 μg/L and by steps of 1 μg/L between 10 and 20μg/L. For any given value of this threshold (say 8 μg/L) the sensitivity was computed as the number of simulated USS-specific urinary exceeding 8 μg/L, divided by 12x200=2400, the simulations corresponding to an atmospheric exposure equal to or exceeding the 20 μg/m 3 TLV. Similarly the specificity was computed as the number of simulated USS-specific urinary below 8 μg/L, divided by 15x200=3000, the simulations corresponding to an atmospheric exposure below the 20 μg/m 3 TLV. Table 3.1 shows the areas under the curve for each of the strategies. As the these areas are between 0.70 and 0.76, there is no great difference between strategies. Strategy D however is virtually identical to the strategy recommended by the ACGIH, so that the fact that we found it to be the best is reassuring.

  in steps of 1 μg/m 3 , plus 35 and 40 μg/m 3 . Thus 12 of the 28 exposure simulated exceeded the 20 μg/m 3 TLV. For each of the 200 simulated subjects and 28 possible exposure levels, five urinary sampling strategies (USS) were compared: (A) Three hour urines collected each day at the shift end, and averaged over four workdays; (B) Three hour urines collected each day after 4 hours of exposure, and averaged over four workdays; (C) First urine of the last day of the workweek; (D) Urine from the last 3 hours of the last shift of the week; (E) Urine from the 3 hours following the end of the last shift of the week.

  urinary cobalt [µg/void]

Figure 2 . 11 :

 211 Figure 2.11 : Scatter plot of observed and predicted urinary data .

  atmospheric cobalt [µg/m3]

Figure 2 . 12 :

 212 Figure 2.12 : Scatter plot of observed and predicted atmospheric data

Figure 2 .Figure 2 . 13 :

 2213 Figure 2.13 : Fifty urinary simulations based on the subject E atmospheric measurements,

Figure 3 .Figure 3 . 1 :

 331 Figure 3.1 shows the graphical comparison of these ROC curves. Five urinary

TABLE I .

 I Summary statistics of the atmospheric and urinary cobalt concentration data in µg/g crea and in µg/L.

			Atmospheric cobalt measurements					Urinary cobalt measurements
	Subject	Number of days	Number of samples	GM (A) µg/m 3	GSD (B)	Min µg/m 3	Max µg/m 3	Number of days	Number of samples	GM (A) µg/g crea	GSD (B)	Min µg/g crea	Max µg/g crea	GM (A) µg/L	GSD (B)	Min µg/L	Max µg/L
	A	4	8	4.9	1.4	4	9	7	35	1.6	2.1	0.5	21.9	2.3	2.3	0.5	22.7
	B	4	8	9.4	1.6	4	18	7	34	3.6	1.7	0.6	8.4	5.7	2.2	1.2	20.7
	C	4	8	10.7	1.4	7	16	7	34	2.9	1.9	1.0	7.8	3.9	2.2	0.8	16.4
	D	4	8	6.2	1.8	4	17	7	33	2.5	1.8	0.8	6.2	2.7	2.2	0.2	10.1
	E	4	8	61.6	2.8	22	276	7	34	12.3	2.3	1.9	87.9	14.3	2.7	2.0 197.8
	F	4	8	26.3	1.2	20	32	7	36	6.0	1.6	2.5	13.1	8.1	2.2	1.4	30.2
	G	4	8	25.6	2.0	10	74	7	35	4.2	2.1	1.2	16.4	4.9	2.2	1.2	38.4
	1	5	10	144.2 2.0	66	449	7	34	16.4	1.9	5.4	65.1	14.5	1.9	1.8	47.7
	2	4	4	26.9	3.4	6	72	4	25	11.5	2.2	1.9	55.1	15.5	2.3	2.0	71.0
	3	6	6	41.8	2.9	17	228	7	34	8.6	2.0	3.0	69.7	5.7	2.4	0.3	34.8
	4	6	6	74.5	2.8	21	205	7	34	11.4	1.6	5.2	28.6	15.0	1.5	7.1	43.5
	5	5	5	14.8	2.9	6	72	4	25	4.5	1.5	2.1	12.7	6.6	1.7	2.0	15.8
	6	5	5	19.8	1.3	14	25	7	29	8.7	1.6	2.9	20.2	9.6	1.5	3.7	22.1
	7	5	5	26.8	1.8	18	77	4	26	7.5	1.7	2.2	13.8	7.2	1.5	3.2	16.6
	8	5	5	23.3	1.6	12	35	4	25	7.2	1.7	2.6	15.9	10.0	1.3	5.6	19.1
	9	6	6	5.9	2.0	3	19	7	34	3.0	1.7	1.1	8.0	2.6	1.9	0.3	6.9
	(A) GM = Geometric Mean
	(B) GSD= Geometric Standard Deviation

TABLE II .

 II Summary statistics of the posterior distributions for population and individual parameters.

	Parameter (unit)	Population geometric mean µ (95% CI)	Between subject geometric SD Σ (95% CI)	Minimum of mean subject-specific value (subject)	Maximum of mean subject-specific value (subject)
	K in (L/h) (A)	60.4 (27.3-97.3)	2.04 (1.47-4.57)	29.9 (1)	123 (D)
	K e (1/h) (B)	0.068 (0.056-0.083)	1.16 (1.02-1.54)	0.059 (F)	0.076 (9)
	T e1/2 (h) (C)	10.19 (8.35-12.37)	1.16 (1.02-1.54)	11.75	9.12
	K r (1/h) (D)	0.078 (0.040-0.25)	1.93 (1.15-4.89)	0.054 (G)	0.171 (8)
	T r1/2 (h) (C)	8.88	(2.77-17.32)	1.93 (1.15-4.89)	12.8	4.0
	K s (1/h) (D)	0.034 (0.013-0.10)	2.22 (1.17-4.98)	0.015 (3)	0.117 (1)
	T s1/2 (h) (C)	20.38 (6.93-57.76)	2.22 (1.17-4.98)	46.2	5.92
	(A) K in : intake coefficient (B) K e : elimination coefficient (C) T x1/2 = half time = ln(2)/K x (D) K r and K s : transfer coefficients between central and peripheral compartment, (see online appendix 2 and text).

Table 3 .

 3 1 : Areas under the ROC curves corresponding to USS

	Urinary Sampling	Areas under
	Strategies (USS)	the ROC curves
	(D)	0.76
	(E)	0.75
	(A)	0.73
	(C)	0.71
	(B)	0.70

Toxicokinetic model

For each subject, the time course of cobalt excretion (in µg) in urine noted Q u , was modelled in a deterministic two-compartment model with central and peripheral compartments. Figure 1.1 presents that model graphically: arrows represent the flows between compartments or the external environment. K in (in L/h) conditions the quantity of cobalt inhaled (C in in µg/L) entering the organism. K e , K r and K s are transfer coefficients (in 1/h) and describe the flow between the compartments (K r , K s ) or excretion (K e ). The corresponding half-times can be computed as ln(2)/K i were i is either r, s or e. The following differential equation system describes the temporal evolution of the quantities of cobalt quantities Q c , Q p and Q u in the central compartment, peripheral compartment and urine, respectively: The input parameter K in , the transfer rate constants K r and K s , and the excretion rate constant K e have a specific value for each subject.

Priors, model checking, cross validation

Priors and posteriors

The Bayesian framework in which our model was calibrated requires to specify priors for its parameters.

Model checking

The fit of the model to the data was examined by plotting the observed vs predicted urinary and atmospheric measurements.

The resulting Figures 2.11 and 2.12 are shown below .


Isabelle Contrucci
email: isabelle.contrucci@ineris.fr

Emmanuelle Klein
email: emmanuelle.klein@ineris.fr

Pascal Bigarre
email: pascal.bigarre@ineris.fr.

Armand Lizeur

Anthony Lomax
email: anthony@alomax.net.

Mustapha Bennani

P Bigarré

Management of post-mining large-scale ground failures : blast swarms field

Keywords: mine collapse, risk management, microseismicity, early warning system, 3D location, velocity model

est destinée au dépôt

Introduction

In the Lorraine area of Eastern France, decades of iron-ore mining from 1850 to 1997 have left vast underground cavities beneath or in the vicinity of urban areas. At present, these residual voids are estimated at 500 million m 3 and represent 40,000 km of underground galleries. This now poses a societal and environmental challenge: major collapses took place in the 1990s in the southern part of the Lorraine iron-ore basin in the cities of Auboué, Moutiers and Roncourt. These events occurred a few months after the mine closure and the progressive rise of the water level in the underground working caused by the halt of the de-watering system [START_REF] Didier | The French Experience of Post Mining Management[END_REF].

These large scale ground failure events prompted a request from both the government and local authorities for a management strategy to be set up to prevent and control post-mining risks. A methodology was thus established in order to assess hazard zones and rank them according to their vulnerability due to human surface infrastructure and activity. The high risk zones are secured either by reducing the hazard or by using in situ monitoring. As it is rarely possible to deal with the hazard itself, public safety is often ensured through real-time microseismic monitoring systems that are installed in the zones in question.

For that purpose, INERIS deployed an innovative microseismic monitoring platform that was tested and validated during the Terres Rouges experiment in 1997 [START_REF] Senfaute | Qualification of the microseismic monitoring technique applied to the risk of collapse or iron mines[END_REF][START_REF] Couffin | Permanent real time microseismic monitoring of abandoned mines for public safety[END_REF][START_REF] Bennani | Coverage formations at the level of areas in which there is a sudden collapse hazard[END_REF]. This platform is designed to detect rock mass fracturing that first affect the old mine workings before reaching the overburden and finally the surface. More than thirty real-time microseismic monitoring networks have been installed since 1998 in the Lorraine iron-ore basin by the CENARIS 1 at INERIS. These networks have been designed following a methodology based on the geological context, the dimensions of the zone being monitored, the hazards involved and the technical installation and maintenance constraints due to urbanization. The most vulnerable buildings are instrumented with geotechnical sensors linked to the microseismic system to ensure a smart triggering scheme [START_REF] Contrucci | Early-warning microseismic systems applied to the management of post-mining large-scale ground failures: calibration by a geophysical field experimentof blast swarms[END_REF].

Although no major collapse has occurred in any of these areas yet, small microseismic events are often recorded. These events, with low amplitudes (~10 -6 mm/s) and negative magnitudes raise real quantitative analysis problems, especially in terms of localisation and spatial-temporal breakdown which are essential in understanding instability mechanisms [START_REF] Driad-Lebeau | Seismic and geotechnical investigations following a rockburst in a complex French mining district[END_REF]. Although publications exist on the microseismic monitoring of working mines (e.g. [START_REF] Senfaute | Spatial distribution of mining tremors and the relationship to rockburst hazard[END_REF][START_REF] Driad-Lebeau | Seismic and geotechnical investigations following a rockburst in a complex French mining district[END_REF][START_REF] Li | A review of mining-induced seismicity in China[END_REF] or on the long-term analysis of microseismicity in abandoned mines [START_REF] Miller | Microseismic techniques for monitoring incipient hazardous collapse conditions above abandoned mines[END_REF][START_REF] Ogasawara | Microseismicity induced by heavy rainfall around flooded vertical ore veins[END_REF], none of them describe the applied methods, tools or data management strategy. This lack of methods, combined with the flooding of the Nord iron-ore basin scheduled to start at the end of 2005, motivated a large-scale field experiment to help to calibrate data processing parameters and to estimate microseismic source parameters. This experiment was performed in six instrumented zones where mine workings were still accessible. It consisted of numerous blasts sequences to ensure the high sensitivity of the monitoring devices and to calibrate some of the fundamental numeric data processing procedures using accurately-controlled input data. This article describes the sequence of the 2005 calibration experiment and the methods adopted to confirm the high sensitivity of sensors based on blast data. Details will then be given on the 1 French National Monitoring Centre for Ground and Underground Risks created by INERIS calibration procedure for the 3D localisation tool and the determination of its main input parameters, especially the velocity model. An empirical law to calculate the source energy from the hypocentral distance, in the geological context of the Lorraine iron-ore basin, is also estimated. Finally, the transfer of these results to the operational point of view will be discussed.

Geological context of the Lorraine iron-ore basin

The Lorraine iron-ore basin, located on the eastern boundary of the Paris basin, extends for approximately one hundred kilometres from North to South and twenty to thirty kilometres from East to West. The Pont-à-Mousson anticline splits the region into two distinctive zones: the Briey-Longwy-Thionville basin to the North and the Nancy basin to the South. The Lorraine iron-ore basin formations are marine in origin and were subject to very little deformation after their deposit.

The basin is nevertheless cut by major faults of several kilometres long oriented NE-SW. The ironore series and its cover show a dip of several degrees towards the South-West, except near the major fractures where the local dip can be as much as 10°, as for example at the Audun-le-Tiche fault.

The iron-ore series belongs to the Toarcian and Aalenian stages. The marly-carbonated cover belongs to the Bajocian. Overall, the various geological formations found in the basin are as follows [START_REF] Bennani | Coverage formations at the level of areas in which there is a sudden collapse hazard[END_REF] The iron-ore formation can be considered as alternating between marls and ores. It is between 0 and 65 m thick, 40 m on average [START_REF] Montagne | Les mines de fer de Lorraine[END_REF]. Nine ore-bearing layers have been counted in the entire basin [START_REF] Maubeuge | Geological observations to the east of the Paris basin[END_REF], named using colours. Usually two to three layers were exploited with the rooms and pillars operating method below urban areas, with pillar extraction outside urban areas. A total of 3.1 billion tonnes of iron ore were extracted, i.e. about 1.2 billion m 3 over a surface area of 1700 km 2 .

Field experiment design

Experimental sites

The field experiment took place between October 2005 and May 2006 in the cities of Audun-le-Tiche, Fontoy, Tressange, Nondkeil, Ottange and Moutiers (Figure 1). These municipalities were already equipped with permanent microseismic monitoring networks; the underground workings were accessible during the experiments on all sites except Moutiers. Each monitoring network includes one or more microseismic stations. Each station consists of three probes equipped with miniature broadband geophones and cemented into boreholes as follows: one 1D probe on the surface, one 1D probe about 15 m below the surface and one 3D probe about 50 m below the surface. The orientation of the 3D probe was measured at the installation into the borehole. The local coordinate system used is the Lambert 1. The stations are connected to an acquisition unit, which automatically detects, records and transfers data to the monitoring central site of INERIS at Nancy, in quasi real-time via, a secure high bandwidth link.

Experimental procedure

Apart from the Moutiers site, where the blasts were made from a vertical borehole, the experiments involved small dynamite blasts in mine pillars at depths of between 50 and 250 m depending on the site. TITADYNE AG 30 dynamite was the explosive used, packed in a cartridge 50 mm in diameter placed in 76 mm-diameter boreholes 5 m deep. The explosive charges were chosen to produce signals with usable signal-to-noise ratio yet avoiding sensor saturation. Thus, the sources were dynamite charges between 0.5 to 12 kg, depending on the experimental site and its configuration.

The blast positions were chosen considering both the best "blast-to-station" distance coverage and the best angular coverage when mine accesses so permitted. Various blast configurations were also used on each site to test the resolution of the numerical processing tools, including 3D location and source parameters calculations. The first combination, a so-called "single blast", involved one single blast in a borehole (Figure 3-a). The second combination, a so-called "orientation blast", comprised two blasts of equivalent charge in the same pillar placed in perpendicular boreholes (Figure 3-b).

Lastly, a combination of orientation and single blasts -"multiplet blast"involved four blasts using different charges in the same pillar, including one in a borehole perpendicular to the other three (Figure 3-c). When possible, fixed X and Y coordinate blasts were also reproduced in other exploited levels to test the robustness of the localisation at depth.

The "multiplet blast" configuration was performed to record "multiplets" a term which refers, in classic seismology, to a group of seismic events showing the following characteristics [START_REF] Geller | Four similar earthquakes in central California[END_REF][START_REF] Poupinet | Monitoring velocity variations in the crust using earthquake doubletsan application to the calaveras fault, California[END_REF][START_REF] Lees | Multiplet analysis at Coso geothermal[END_REF][START_REF] Slunga | Absolute and relative locations of similar events with applisation to microearthquakes in southern Iceland[END_REF][START_REF] Moriya | Analysis of fracture propagation behavior using hydraulically induced acoustic emissions in the Bernburg salt mine, Germany[END_REF][START_REF] Gibowicz | Seismic doublets and multiplets at Polish coal and copper mines[END_REF]: similar wave form; different time origin; considered as the result of the relaxation of constrain of a single fracture or fault, similar source mechanisms; close location; different magnitude. The recorded data might indeed be used to test the capabilities of "relative" location algorithms compared to "absolute" and classical seismic location algorithms. The relative approach is supposed to improve the location of one order of magnitude compared to absolute location [START_REF] Rubin | Streaks of microearthquakes along creeping faults[END_REF][START_REF] Waldhauser | A double-difference earthquake location algorithm: Method and application to the Northern Hayward Fault[END_REF][START_REF] Schaff | Highresolution image of Calaveras Fault seismicity[END_REF].

Seventy blasts were performed in total, producing over 1270 seismograms (Table 1). The recorded seismograms show a frequency spectrum of between 30 and 500 Hz for hypocentral distances of less than 300 m (Figure 4). This frequency range was found to "fit" correctly the frequency range observed from a previous experiment during which caving induced events were recorded [START_REF] Senfaute | Qualification of the microseismic monitoring technique applied to the risk of collapse or iron mines[END_REF].

Site

Data processing: example of the Tressange site

Microseismic data processing started with a systematic analysis of all the seismograms in terms of amplitude, signal-to-noise ratio and sensor energy to assess data quality. All the recorded 3D seismograms were then analysed for polarisation to check consistency between measured incident angles and expected values at the 3D probes, thereby ensuring the correct azimuthal orientation of these probes ( §4.2). The geological formation velocities for the overburden were then estimated from the inversion of both P-wave arrival times and measured polarisation angles, taking into account the multilayered geology of the site. The blasts were re-localised based on the velocity model thus calculated ( § 4.3.3). Lastly, an empirical law for calculating source energy was determined based on the sensor distance and energy ( §4.4). This approach is illustrated below with the Tressange check blast data (Figure 2), except for the source energy calculation where all the sites are considered.

4.1.Data quality

A systematic analysis of all seismograms per site and microseismic station in terms of amplitude and energy was undertaken in order to assess data quality. This analysis demonstrated the excellent sensitivity of the monitoring systems, able to detect 1 kg blasts at a distance of up to 300 m. For analysis, only signals with signal to noise ratio higher than 10 were considered. Note that in such a "near-field" experiment, where sources generate mainly compressional stress waves, only P-waves are observed; the S-wave energy is not observed on the recorded data. This limitation will be discussed here below.

4.2.Polarisation analysis

Analysing polarisation by wave rotation is used to calculate incident angles (azimuth and dip) of the ray path at the three-component probe from a 3D seismogram. The estimated direction of the incident ray path provides important information that can be used along with the measured P and S arrival times to constrain efficiently the solution for source location. This approach, which assumes prior knowledge of the orientation of the 3D probes, is of most importance when localisation relies on only a few microseismic stations [START_REF] Magotra | Seismic event detection and source location using single-station (three-component) data[END_REF][START_REF] Abdul-Wahed | Source location estimation using single station threecomponent seismic data[END_REF][START_REF] Volker | Automated seismic event location for hydrocarbon reservoir[END_REF]. In the current study, the azimuths of the 3D probes were measured with a releasable compass during grouting into the vertical boreholes. As the blast positions are known, the orientation measurements can be checked easily by examining the expected and measured directions by a polarisation analysis of each 3D seismogram. The expected direction is calculated by assuming an isotropic homogeneous medium between the blast and the 3D probe, i.e. assuming a straight ray path between these two points. This hypothesis is reasonable regarding the short "blast-to-probe" distance and the relative homogeneity of the encountered geology. Note that in such sub-horizontal media, the calculated azimuth does not depend on the velocity model. On the contrary, the measured dip of the ray path may depend significantly on the gradient velocity of the stratified overburden when contrasting velocities between layers are considered.

On the Tressange site, the expected azimuths at the 3D probes at the Nationale and Liberté stations match the measured azimuths; the mean error, assessed using quality factors such as rectilinarity and planarity factors [START_REF] Samson | The spectral matrix, eigenvalues, and principal components in the analysis of multichannel geophysical data[END_REF], is ~5°. However, for the Gerbault and Jardins probes, a systematic mean misfit of 50° and -200° respectively is noted between the expected and calculated azimuths (Figure 5). These systematic misfits are caused uncontrolled rotation of 3D probes during grout injection while the orientation device is released and pulled out. These misfits have been corrected at each station, thus providing true probe orientation. Note that blast 08, located underneath the Liberté station, has a major azimuthal misfit (Figure 5) because of the inconsistent solution for ray path with a dip close to 90°. Note also that a bias between expected and measured dip values is observed for all stations: the bias is less than 5° for the Liberté, Jardins and Nationale stations (Figure 6), it is ~ 9° for the Gerbault station (Figure 6). It is most likely due to straight ray paths assumption.

3D localisation of the blasts

This stage consists in performing blast re-localisation to test the performance of the localisation algorithm thanks to the known blast positions ( §4.3.4). Firstly, it is necessary to build up a velocity model based on the site's geological structure ( §4.3.1). P-wave velocities for each layer are calculated taking account all available input data, i.e. positions of the blasts, direct P-wave arrival times and dips measured at the 3D probes (see §4.2). Note that S wave arrival times were not available and measured azimuths are of no interest in the velocity calculation. For this purpose, a velocity model optimization program has been developed (4.3.3).

Geology of the Tressange site

The Tressange site is located on the Eastern part of a fault system extending from Mont-Bonvillers to Ottange, oriented SW-NE to SSW-NNE. Note the presence of the Ottange fault, orientated SSW-NNE, at Tressange (Figure 2). The fault throw is several metres and the East block is the collapsed block with layers dipping regionally in the order of 3% westwards or WSW [START_REF] Lopes | Etude de l'environnement géologique du site de Tressange[END_REF]. As shown in Figure 2, the Nationale and Jardins stations are situated in the West fault block whereas the other stations are located in the East block. The Ferdinand mine shaft (Table 2), which is located to the west of the Gerbault station, in the East block, shows that the geological structure comprises six main facies (Table 2). The geological structure adopted for Tressange is based on these different characteristics.

Coordinates

Construction of an inversion type velocity model

SYTMISvel software was specially developed for automatically adjusting a velocity model to a determined geological structure comprised of inclined parallel layers with known thicknesses [START_REF] Contrucci | Early-warning microseismic systems applied to the management of post-mining large-scale ground failures: calibration by a geophysical field experimentof blast swarms[END_REF]. This model is build up by inversion, using the arrival times of the P and/or S waves and the polarisation angles. Velocities are calculated with or without knowledge of the absolute initial time (T0) of the recorded seismograms.

SYTMISvel program explores numerous velocity models, which have been generated randomly, according to a Monte-Carlo type algorithm [START_REF] Lomax | Identifying sets of acceptable solutions to non-linear, geophysical inverse problems which have complicated misfit functions[END_REF]: the applicability of these models is assessed according to the differences between calculated and observed data with respect to their uncertainty. The models achieving the most effective minimisation of the difference between calculated and observed data are selected, if the velocities calculated in the individual layers are consistent with the local geology.

The differences between calculated and observed data are minimised by using either the least-square misfit function (L2-norm) or the Equal Differential Time misfit function (EDT-norm; [START_REF] Pinsky | A comparative study of robust algorithms for rapid, automatic earthquake location[END_REF][START_REF] Lomax | A reanalysis of the hypocentral location and related observations for the great 1906 California earthquake[END_REF].

For the L2-norm, the function to be minimised is expressed as:

j i obs j j j obs i i i LSL x calc x obs x Tcalc x Tobs k pdf 2 2 2 2 2 ) ( ) ( ) ( ) ( exp
with: x, the position of the source; Tobs i and Tcalc i , the arrival times observed and calculated from observation; obs i and calc i are the polarisation angles observed and calculated; i and j are the uncertainties; k a normalisation constant.

For the EDT-norm , the function to be minimised is expressed as: This minimization process uses a global approach since the travel times are non-linearly related to the velocities in the individual layers. The program input comprises mainly the following information: 1) the method used to randomly generate sets of velocity models either the traditional Monte-Carlo method which defines randomly the entire velocity model or the "Metropolis" Monte-Carlo method in which each model tested in a Metropolis loop, is close to the last model selected (a random walk similar to the simulated annealing method); 2) the geological structure and the observed data, i.e., the number of geological layers and the velocity range for P and S-waves; the thicknesses of the layers and their possible dips; the positions of the microseismic stations and blast points; the arrival times of the P and/or S waves for each blast and each of the 1D and 3D probes, as well as associated picking errors and the incidence angles observed for each blast at each 3D probe and associated errors. The program output comprises principally the difference between calculated and observed arrival times and incidence angles for each blast and each of the considered probes and the probability of occurrence of each velocity model.

) ) ( ) ( ) ( ) ( exp ) ( ) ( ) ( ) ( exp ( 1 , 2 ,

Velocity model for the Tressange site

The thicknesses of the model's layers are based on the geological structure observed on the log of the Ferdinand Shaft (Table 2). Several tests were carried out in order to determine the best velocity model with the layer thicknesses indicated in Table 2. The minimisation between calculated and observed values was achieved by using the L2-norm. Regardless of the configuration of the data used for the inversion, significant misfits were indeed obtained when the EDT-norm was employed. This norm does not seem to be adapted to the geometry of the problem and the dimensions of the geological structure studied: the double difference of the arrival times and angles affects the information contained in our data, and thus generates aberrant velocity values.

The best velocity model is determined in two main steps. The first steps consists in inverting all of the data, i.e. all of P-wave arrival times and incidence angles of all blasts. It allows identifying the data with large misfit between calculated and observed values (Figure 7-a-b-c). Figure 7-a shows that time misfits higher than 0.04s are observed for blast 1 on Gerbault station, and blasts 1, 4 and 14 on Liberté station. The high misfit observed for blast 1 on Gerbault station is probably related to a local heterogeneity in the vicinity of the blast point. This hypothesis seems reasonable since the hypocentral distance is, in that case, short. On the contrary, the high time misfits observed for blasts 1, 4 and 14 on Liberté station are probably due to the large hypocentral distances. These high distances can induce significant signal attenuation, leading to high picking errors. These misfits and errors can also be emphasized by anisotropic wave propagation. Figure 7-b shows three significant azimuth misfits. For blast 8, at Liberté station, the misfit is due to the fact that the blast is located directly below the Liberté station. Indeed, there exists an infinity of azimuth solutions for a dip close to 90°. For blasts 4 and 14, the azimuth misfits observed respectively at Gerbault and Nationale stations, are difficult to explain since the hypocentral distances are short in both cases. These misfits are however lower than 10%. They can be neglected since they do not depend on the velocity model.

Regarding the dip (Figure 7-c), the greatest misfits are observed for large hypocentral distances, relatively large error in dip for near-horizontal rays in the assumed, constant velocity layered model.

As a general rule, the highest misfits, in terms of both polarisation angles and P-arrival times, are observed for blasts performed on the network border. This is probably due to signal attenuation with distance. Also, the boundary ray paths cross geological sequences that are not intersected by other ray paths: the information collected from the boundary ray paths is not crosschecked thus creating high misfits. In other words, this first stage demonstrates that the velocity model has to be calculated with high signal to noise input data. Thus, the second step in the velocity model calculation was carried out by eliminating the "aberrant" errors or high misfits (Figure 8) associated with blasts 1, 4 and 14.

Step The step 2 model significantly reduces misfits which exist between calculated and measured values for P wave arrival times as well as for the polarisation angles of the incident ray (Figure 8). Thus, time misfits vary in a maximum range of 0.02 s instead of 0.07 s after the first step inversion.

The velocities obtained by inversion, for each layer, are indicated in Table 3. The results show that removing aberrant input data considerably reduces the standard deviations, except for Jaumont Limestones which layer is clearly not constrained enough by the input data. Figure 9, which presents the distribution of velocities for each layer, gives additional details on these results. The velocity peak of the Polypiers Limestones layer shows that this layer is very well constrained. For the other layers, the peaks are less well-defined (Figure 9), and thus associated errors are greater (Table 3).

These calculated velocities can also be compared with velocity measurements made in laboratory by [START_REF] Homand | Survey of the coverage formations of areas in which there is a non-negligible risk of sudden collapse, LAEGO, LG[END_REF]. These measurements give an accurate idea of expected velocities in the formations in question, and, the contrast between layers. They show that P-wave velocities vary between 2710 and 4215 m/s for the Polypiers Limestones, between 3900 and 5375 m/s for the Ottange/Haut-Pont Limestones, between 1965 and 2680 m/s for the Charennes Marls and between 3270 and 3450 m/s for the iron-ore formation. These values are consistent with those obtained in this study. Note that the step 2 velocity model shall be selected for re-locating the blasts.

3D localisation algorithm

The localisation module implemented in SYTMISauto software used in this study, is based on the combination of the microseismic wave arrival times as well as the polarisation angles, in order to determine the hypocenter with the maximum likelihood. Indeed, the integration of polarisation angles enables to locate an event with few probes, i.e. one 3D probe and one 1D probe if only P waves are detected or a single 3D probe if P and S waves are recorded [START_REF] Magotra | Seismic event detection and source location using single-station (three-component) data[END_REF][START_REF] Abdul-Wahed | Source location estimation using single station threecomponent seismic data[END_REF][START_REF] Volker | Automated seismic event location for hydrocarbon reservoir[END_REF]. The implemented localisation algorithm is based on a probabilistic approach to solve the inverse problem [START_REF] Tarantola | Generalized Nonlinear Inverse Problems Solved using the Least Squares Criterion[END_REF]. It consists in maximising the probability density function (pdf) of the hypocenter at a given point using EDT or L2 norm. This is done by minimizing the misfit between measured and calculated values, i.e. between observed and calculated arrival times, as well as between observed and calculated polarisation angles.

The hypocenter with the maximum likelihood is determined by using the Oct-Tree non-linear method [START_REF] Lomax | Fast, probabilistic earthquake location in 3D models using oct-tree importance sampling[END_REF], based on a successive division of space into cells depending on the value of the probability calculated for each cell. This approach is used to completely solve the inverse problem and thus provides a representation of the overall pdf of the localisation. The most probable hypocentre corresponds to the pdf maximum. The drawback of this method is that the calculation can take a long time because the entire solution space is explored for each iteration. To limit time calculations, this method can be combined with the linear C.H.E.A.P method [START_REF] Tarantola | Generalized Nonlinear Inverse Problems Solved using the Least Squares Criterion[END_REF] based on the determination of a local maximum by gradient calculation, that is to say of the maximum pdf in the relevant cell. For each event, the pdf, as well as the corresponding reliability ellipse at 68 %, is calculated.

Localisation parameters

Calibration blasts are re-located on the Tressange site with the global Oct-Tree method to avoid the restrictions of the linear C.H.E.A.P method described above. The grid size was set as 10 m minimum and 50 m maximum. In order to optimise the calculation time, the grid space was limited in a 1500 m * 1500 m * 500 m cubic space following X Y Z directions. As for the velocity model calculation, localisation processing is run based on the L2-norm, since the EDT-norm did not seem to be adapted to the dimensions of the studied geological structure. Although it seems appropriate for locating regional or worldwide earthquakes by accounting for actual picking errors [START_REF] Pinsky | A comparative study of robust algorithms for rapid, automatic earthquake location[END_REF][START_REF] Lomax | A reanalysis of the hypocentral location and related observations for the great 1906 California earthquake[END_REF], the EDT-norm generates aberrant localisation solutions in the current study.

The velocity model employed is from SYTMISvel software ( § 4.3.3). To partially take into account uncertainties related to the velocity model determined previously by inversion, a constant error with value of 0.002 s is introduced on all travel times, which corresponds to an error of 70 m/s for a stress wave travelling at 3500 m/s over a distance of 350 m. Errors on the incidence angles are assessed at 10° which seems reasonable according to Figure 8; manual picking error is set at 0.005 s.

3D localisation results

Several re-localisation tests were carried out on the 13 blasts related to the site of Tressange. Firstly, only P-wave arrival times were used; secondly, the polarisation angles (azimuths and dips) were added. The results of these tests are summarised in Table 4 and Figure 10. When only P wave arrival times are considered, the mean misfit between the true and calculated positions of the blasts is 260 m with significant errors in the Z direction, i.e. depth (210 m on average). When P-wave arrival times and polarisation angles are taken into account in the localisation, the mean misfit between actual positions and calculated positions is reduced to 70 m (Table 4). There is a misfit of 48 m in the X -Y plane and of 47 m in the Z direction (Table 4). The localisation quality is improved in nearly every case, except for blast 08 located directly below the Liberty station. In this specific case, the azimuth measurement has little meaning as it is determined with a very significant error, which affects localisation quality. On the other hand, blasts 09, 10, 11 and 12 located in the centre of the network are located correctly with an error of about 55 m on average. For blasts 01, 04, 13 and 14, located outside the network, signal to noise ratio of the signals recorded on the distant stations is too low for clearly for accurate picking of the first arrivals. The localisation of these blasts, with only two or three stations, is greatly improved with the polarisation angles (Table 4 and Figure 10).

More generally, when only the arrival time is considered, the pdf reveals a North West -South East orientation, related to the geometry of the network. The pdf is far more constrained when polarisation angles are taken into account (Table 4). However, both the input data and the precision of the 3D localisation algorithm remains insufficient for the data recorded with multiplet-type configuration (blasts 9, 10, 11 and 12). It would certainly be useful to use relative localisation methods [START_REF] Poupinet | Monitoring velocity variations in the crust using earthquake doubletsan application to the calaveras fault, California[END_REF][START_REF] Slunga | Absolute and relative locations of similar events with applisation to microearthquakes in southern Iceland[END_REF][START_REF] Abdul-Wahed | Mining-induced seismicity: Seismic measurement using multiplet approach and numerical modelling[END_REF] to gain precision while re-localizing this kind of data. 195560 195504 195456 195449 195268 195129 195126 195122 195115 194914 195409 195409 195245 195022 195034 195022 195022 194799 195514 195491 195426 195426 195245 195092 195098 195092 195092 194864 Table 4 : Summary of the true and re-calculated positions of microseismic blasts, hypocentral misfits in the X-Y plane and depth (Z direction), and associated average values.

4.4.Seismic energy at the source

The seismic energy at the source is a key parameter for characterising the intensity of a seismic rupture. It is however necessary to assess the geometrical and anelastic attenuations of the rock mass and the response of the instrumental channel [START_REF] Urbancic | Space-time variations in source parameters of mining-induced seismic events with M < 0[END_REF][START_REF] Boatwright | Regional Estimates of Radiated Seismic Energy[END_REF][START_REF] Yamada | Stress drops and radiated seismic energies of microearthquakes in a South African gold mine[END_REF]) in order to obtain a correct estimate of this energy. Although it is easy to correct the recorded signal of the geometrical attenuation, it is more difficult to perform other corrections without making numerous assumptions. In conventional seismology, the calculation procedures used for characterising the source often involve a double-couple mechanism since fault shearing mechanisms usually generate earthquakes [START_REF] Brune | Tectonic stress and the spectra of seismic shear waves from earthquakes[END_REF]Madriaga, 1976). Based on this approach, the magnitudes of the earthquakes can be assessed, thus allowing the determination of the source energy [START_REF] Aki | Quantitative Seismology[END_REF]. However, these techniques cannot be applied directly for assessing the source energy of a micro-seismic event occurring in a mine. Source mechanisms, different from the double-couple, are expected, for instance, traction, implosion or explosion [START_REF] Gibowicz | Source Parameters of Seismic Events at the Underground Research Laboratory in Manitoba, Canada: Scaling Relations for Events with Moment Magnitude Smaller than -2[END_REF][START_REF] Mcgarr | An implosive component in the seismic moment tensor of a mining-induced tremor[END_REF][START_REF] Trifu | The use of uniaxial recordings in moment tensor inversions for induced seismic sources[END_REF][START_REF] Finck | Advances in Moment Tensor Inversion for Civil Engineering[END_REF][START_REF] Sileny | Seismic Moment Tensor Resolution on a Local Scale: Simulated Rockburst and Mine-induced Seismic Events in the Kopanang Gold Mine, South Africa[END_REF].

Analyses are currently being conducted to determine and quantify the different corrections to be taken into account in the calculation of the source energy in the context of this study. Therefore, this paper only includes a description of the work concerning the empirical estimation of the source energy from calibration blasts. This approach appears to be a good first order solution since the blasting charges and their positions are known. An empirical relationship between the seismic energy recorded by the sensors and the seismic source energy according to the distance has been defined as follows [START_REF] Tastet | Large-scale field experiment to calibrate microseismic source parameters applied to real-time monitoring of post-mining instabilities[END_REF]:

D Esrc K Ecpt (1)
with Esrc: source energy linked to the explosive charge (J), see equation ( 2); Ecpt: energy recorded by the sensor and calculated from seismograms (J); D: hypocentral distance based on straight ray path assumption (m); K: constant (J 1-m ). Note that in such a formulation (1), unknowns are K, and where quantifies both geometrical and anelastic attenuations and dissipation and nonlinear source size effects.

The seismic source energy depending on the blasting charge can be expressed as follows:

r Etotal Q Esrc
(2)

with: E src , seismic source energy (J); Q, source blasting charge (kg) ; E total , total energy released by a blasting charge of 1kg depending on explosive used; r, ratio set empirically at 0.6.

The value of the total energy released (E total ) by a blasting charge of 1 kg depends on the blasting agent being used. In this case, it is TITADYN AG 30 which is packaged in cartridges of 50 mm and 38 cm in length. As these calibration blasts have been carried out in boreholes (confined space) i.e. in a "blocked blast" configuration, the maximum detonation velocity of 6,000 m/s has most probably been reached as well as the maximum energy of 4.2.10 6 J for 1 kg of blasting charge. This energy includes two terms: the first represents the shock energy and the second the gas energy. Only the shock energy can significantly impact the wave transmission in the rock and contribute to the seismic energy radiated into the rock mass with an insignificant participation of detonation gas. Thus, coefficient r with a value of 0.6 has been applied to the total energy in order to account for this phenomenon. Equation 1 leads to the following expression:

1 * D K Ecpt Esrc (3)
where Esrc*: estimated source energy (J); K, and coefficients are determined, from equation 1, for each experimental site. First, is determined for a given series of blasts of constant source energy Esrc, defined by the blast charge. (considered here as a constant with the blast charge) and power law .

Secondly, once has been determined, K and coefficients are determined to fit all data for the studied site. Eventually, the robustness of the empirical law ( 2) is evaluated from the standard misfit between estimated Esrc* and the true values of Esrc. The results give a mean error of ~50%, which is very acceptable when one considers all cumulated assumptions. The empirical laws obtained for each site are presented in Table 5. Table 5 : Source energy estimation laws from sensor energy for each experimental site.

Site

For each studied site, the dependency of the seismic source energy on the energy recorded by a sensor reveals a good homogeneity with -1 values close to 0.5 except for the Fontoy site where this value is equal to 0.8 (Table 5). This value illustrates the significant attenuation observed on this site where it was necessary to perform blasts of up to 12 kg. The dependency of the source energy on the hypocentral distance reveals relatively homogeneous exponents between 2.3 and 4. The relative homogeneity observed on sites, which are relatively distant from each other, illustrates the geological homogeneity of the Lorraine iron-basin (Figure 11-a) and contributes to the definition of a source energy estimation law common to all sites. In order to develop this law, the same approach employed for a single site has been adopted (Table 5).

This approach, although empirical, can be used for estimating the source energy of the study sites from the measured sensor energy and the position of the event. However, these relationships are valid only for P-waves as the blasts have not generated other types of waves. The influence on these empirical relationships of energy related to S waves, which could be generated by underground failure, must still be assessed. The empirical approach developed here, may however be used to estimate the seismic source energy of real microseismic events on a relative local scale for fast and accurate classification of different events from a unique swarm.

Summary of results and discussion

The experiment conducted allowed the recording of more than 1 200 high quality seismograms, confirming the ability of the microseismic systems deployed to detect signals caused by blasts of 1kg of explosives at a distance of more than 300 m. The data signal to noise ratio quality, along with the polarisation analysis allowed the orientation of 3D probes to be checked and corrected, when necessary, with an accuracy close to 1°. For the Tressange site, two 3D probes over four probes presented a systematic error. These probes no doubt rotated in the borehole. On the whole, this experiment allowed the orientations of 7 out of the 22 3D probes in question to be corrected and brought about a change in the installation protocol of the probes. These are now installed with a centring system, with extra pressure contacts, limiting any accidental rotation.

As already mentioned, accurate knowledge of the 3D sensors orientation is crucial in the processing and analysis of microseismic data. Together with the arrival times measured at the sensors, through the inversion process, angles allow a realistic velocity model, in relation to the site's geology, to be calculated. However, for the Tressange site it was shown that the calculation of the velocity model depends on good knowledge of the geological structure of the studied site. It depends as well on a precise analysis of parameters to eliminate any abnormal data. The use of an appropriate norm to minimize any differences between the observed and calculated data is also important. In the current study, the L2-norm is far more accurate than the EDT standard, for both velocity model calculation and blast re-location.

Calculations of azimuth and dip of the incident rays at the 3D probes is easy and accurate. Once these angles included in the 3D location algorithm, the accuracy of the absolute location hypocenters increases significantly. For the Tressange site, the average location error is 260 m when only the arrival times are used and this is reduced to 70 m when the polarisation angles are also included. This accuracy is sufficient compared to a monitored area of several hectares. This improvement is valid for almost all the blasts, with the exception of blasts located right below stations where the azimuth error is significant. The accuracy of the 3D absolute location algorithm obtained remains insufficient to take advantage of the data recorded in multiplet configurations.

An empirical law providing a relationship between the energy at the source, the energy measured at the sensor and the hypocentral distance was defined for each experimental site. The good homogeneity of these laws reflects the uniformity of the geology at the scale of the Lorraine iron-ore basin. An empirical law, valid throughout the iron-ore field, was estimated. These empirical energy laws are valid for P-waves only, since the blasts did not generate S waves. The influence of any S waves, which may be generated by natural failures, should thus be further evaluated. At this stage, the ability to estimate the source energy of a real microseismic event and compare it to blasts test gives already a good order of magnitude.

Conclusion

The microseismic calibration blasts carried out in the Lorraine iron-ore field form a reference database consisting of 1200 high-quality seismograms. This unique database allows the validation of the tools developed for microseismic analysis and calibration of the microseismic characteristics of each experimental sites.

Step-by-step analysis of those data enabled us to: (1) calculate P-wave velocity models of each site and characterize some geophysical properties of the main geological strata ;

(2) optimize and validate the 3D localization tool; (3) construct empirical source energy and wave propagation laws for the different but geologically similar experimental sites including comparison tests. In the future, we will perform additional research to further characterize the seismic source mechanism. This work will also take into account the relative location method to take advantage of the various blast configurations implemented in this experiment. The results already achieved and the work in progress, should allow the characterisation of post-mining instabilities and their early warning systems in order to enhance decision making.

  : (1) the Doncourt oolitic Limestones; (2) the Jaumont oolitic Limestones; (3) the upper and lower Polypiers Limestones; (4) the Haut-Pont Limestones; (5) the Ottange Limestones; (6) the Charennes Marls; (7) the iron-ore formation. The geometry, thickness and extent of these series can vary laterally rapidly. Not all these formations are found in all the experimental sites. For example, Doncourt Oolitic Limestones are only found in the west of the iron-ore basin.

  a and Tcalc b , the arrival times for observations a and b; Tcalc a and Tcalc b , the calculated travelling times; obs a and obs b are the polarisation angles observed for observations a and b; calc a and calc b are the polarisation angles calculated; a and b are the respective uncertainties; K, a normalisation constant.

Figure

  Figure 11 illustrates the correct determination of ] [ Esrc K(considered here as a constant with the

Figure 2 :

 2 Figure 2: Mining map of the Tressange site. Grey spots represent the microseismic stations; Dark

Figure 3 :Figure 4 :

 34 Figure 3: configuration of a) a source blast, b) an orientation blast, c) a multiplet blast.

Figure 5 :Figure 6 :Figure 7 :Figure 8 :Figure 9 :Figure 10 :Figure 11

 567891011 Figure 5: Azimuth error versus dip for the 3D microseismic probes of Tressange network.

kg -12.0 kg 1278

	Number of microseimic stations	Number of blasts	"Blast-to-station" distance Min (m) Max (m)	Mass of explosives	Number of seismograms
	Fontoy	2	13	100	900	3.0 kg -12.0 kg	130
	Tressange	4	13	170	1000	2.0 kg -5.0 kg	260
	Audun-le-Tiche	3	11	80	1050	2.0 kg -9.0 kg	165
	Nondkeil	4	10	110	370	0.5 kg -2.0 kg	140
	Ottange	5	18	90	1300	2.0 kg -10.0 kg	468
	Moutiers	4	5	80	480	1.0 kg -5.0 kg	115
	TOTAL 0.5 Table 1: table summarising the number of blasts, blast to station distance, charges and signals 22 70 80 1300
	recorded per site.

Table 2 :

 2 Geological facies found at the Ferdinand shaft on the Tressange site.

		(Lambert 2 extended)
	FERDINAND SHAFT
	X = 864456.9 Y = 2495763.8 and Z = 347.51 m
		Depth of the layer	Layer thickness
	Geological facies:	roof (m)	(m)
	Overburden	-	-
	Jaumont Limestones	0	14
	Polypiers Limestones	14	85
	Haut-Pont/Ottange
	Limestones	99	54
	Charennes Marls	153	28.4
	Iron-ore formation	181.4	-

1: all of the data Formation Jaumont Limestones Polypiers Limestones

				Haut Pont/
				Ottange	Charennes	Iron-Ore
				Limestones	Marls	Formation
	Velocity [m/s]	5 530	4 820	4 790	3 235	3 215
	Standard
	Deviation [m/s]	430	215	1 280	1 435	1 510
			Step 2: without blasts 1, 4 and 14
	Velocity [m/s]	2 340	3 375	3 855	3 190	3 190
	Standard
	Deviation [m/s]	600	165	920	970	305

Table 3 :

 3 best P-wave velocity model obtained by inversion in steps 1 and 2.

Acknowledgements

Special thanks to the French Ministry of Economy, Finance and Industry for the financial support.

We fully acknowledge M. E. Rovinalti (former mine manager from ARBED company), M. A. Marnet (DRIRE Lorraine), M. J-P. Josien (former director from GEODERIS) both for authorizations and helpful contribution to the field experiment.

Our thanks also go to the ARBED and ARCELOR companies for the technical support and human resources used during the experiment. We thank the two anonymous reviewers for their suggestions that significantly helped to improve this paper. We do not forget Kathe Hooks (School of Mines of Nancy) for fully re-read this paper.