Cyrille Balland
email: cyrille.balland@ineris.fr

Jacques Morel
email: jacques.morel@andra.fr

Gilles Armand
email: gilles.armand@andra.fr

Will Pettitt

ULTRASONIC VELOCITY SURVEY IN CALLOVO-OXFORDIAN ARGILLACEOUS ROCK DURING SHAFT EXCAVATION

The excavation of a deep underground structure induces a stress field redistribution that could create an Excavation Damaged Zone (EDZ). The study of EDZ is important in the framework of nuclear waste storage, where such EDZ could constitute a preferential pathway of stored material towards the biosphere. Analysis of ultrasonic wave propagation across the rock media around the excavation is one technique that can be used to characterize such zones. This method has been applied in the Meuse/Haute-Marne Underground Research Laboratory (URL) in France to characterise the initialization and evolution of damage around the URL access shaft. To achieve this, a device comprising 21 ultrasonic sensors was spread out at different distances from the shaft wall, at depths of between -465 m and -475 m. A transversal isotropic model enabled the measured experimental velocities on each ray path during shaft excavation to be integrated. The results highlighted the initialization and extension of the damage as the shaft front proceeded. They also showed a polarisation of the velocity field and an oscillation of the transversal isotropy with a preferential orientation of the stress release and the microcracking. Monitoring of the velocity field also highlighted some significantly different behaviour linked to the complex mechanical behaviour of the argillaceous rock.

INTRODUCTION

Various studies in the Canadian URL of Lac du Bonnet [START_REF] Martin | Observations of brittle failure around a circular test tunnel[END_REF] for granite, Swedish HRL [START_REF] Pettitt | The assessment of Damage Around Critical Engineering Structures Using Induced Seismicity and Ultrasonic Techniques[END_REF] for diorite, Mont Terri [START_REF] Thury | The Mont Terri rock laboratory, a new international research project in a Mesozoic shale formation[END_REF], Tournemire [START_REF] Bonin | Deep geological disposal in argillaceous formations: Studies at the Tournemire test site[END_REF] and Mol [START_REF] Mertens | Characterization of induced discontinuities in the Boom Clay around the underground excavations (URF, Mol, Belgium)[END_REF] for argillaceous rock and Asse [START_REF] Hou | Mechanical behaviour of salt in the excavation disturbed zone around underground facilities[END_REF] for salt rock have shown the presence of a damaged zone around the deep excavations. The characterization of such zones is significant in the safety assessment of nuclear waste underground storage. Such zones could constitute a preferential pathway of radionuclide towards the biosphere. Numerical modelling (mechanical, hydromechanical, etc.) provides a first estimation of the extent and the characteristic of the EDZ. However, the results are strongly dependent on the parameters of the model used and the input data. Consequently, the properties of the traversed rock mass need to be well understood through in situ measurements.

Monitoring of ultrasonic wave propagation has demonstrated its capacity to observe the initialization of damage in a previously unperturbed zone. This method is very efficient in granite [START_REF] Carlson | Acoustic emission and ultrasonic velocity study of excavation-induced microcrack damage at the Underground Research Laboratory[END_REF] or in salt rock because wave damping is very low, but it is difficult to apply to argillaceous rock [START_REF] Forney | Caractérisation par méthodes ultrasoniques de la zone endommagée induite par le creusement d'un tunnel en milieu argileux : Cas d'étude au tunnel du Mont Terri[END_REF]. New technology for optimizing the pulsed energy has been developed in the framework of the European OMNIBUS project [START_REF] Pettitt | Development of the Tools and Interpretation Techniques for Ultrasonic Surveys to Monitor the Rock Barrier around Radioactive Waste Packages in Geological Repositories[END_REF]. This technology has been directly applied for monitoring the digging of the main shaft of the Meuse/Haute-Marne URL.

The principle of the method is based on the strong dependency of the pressure wave (P-wave) velocity in the rock with the effective stresses and the damage. The objective of this study is to determine the spatial extension of the damage around the shaft and to highlight its evolution as the shaft excavation front progressed. The device has been designed on the basis of the estimation of the extension of the EDZ by numerical modelling and the intensity of the damage with some preliminary laboratory tests on rock samples.

An analysis method was developed to interpret the velocity changes of the ultrasonic waves in a finite volume in terms of orientation and amplitude. The results have enabled the natural heterogeneity and anisotropy of the rock and the nature and the extension of the couple damage/stress variations induced by the shaft excavation to be characterized.

EXPERIMENTAL LAYOUT 2.1 SITE AND EXPERIMENT DESCRIPTION

The Meuse/Haute Marne URL is located in the north-eastern part of France, at the boundary between the Meuse and Haute Marne Departments. The research activities carried out by ANDRA (French National Radioactive Waste Management Agency) are dedicated to reversible, deep geologic disposal of high-activity, long-lived radioactive waste [START_REF] Delay | Making of underground scientific experimental programme at the Meuse/Haute-Marne underground research laboratory. North Eastern France[END_REF]. The studied rock formation is a 130 m thick clay formation (argillite) lying at a depth between 420 m and 550 m vertical to the URL. Two 500 m deep shafts give access to the research drifts situated at two levels: a first drift at 445 m depth, accessible from the main shaft, and the main URL level at 490 m depth with several research and technical drifts, accessible from the secondary shaft (Fig. 1). The stress state around the URL has been determined from borehole campaigns from the surface [START_REF] Wileveau | Complete in situ stress determination in an argillite sedimentary formation[END_REF]. Principal horizontal stress directions were found to be oriented along N155°E (major stress) and N65°E (minor stress), with a minor stress of about 12.4 MPa and a stress anisotropy of 1.3 at the drift levels. Vertical stress was found to be about 12.7 MPa.

At 445 m depth, the excavation of the main shaft was stopped for the construction of the first drift. From this drift, a specific mine-by experiment, called REP, was carried out. The goal of this experiment was to study the impact of the shaft excavation on the hydromechanical behaviour of the rock mass [START_REF] Su | Hydromechanical coupling phenomena observed during a shaft sinking experiment in a deep argillaceous rock[END_REF] [START_REF] Souley | Poromechanical behaviour of deep claystone and permeability changes around shaft[END_REF]. Several boreholes were drilled from this drift, dipping towards the shaft, and equipped with deformation and pore pressure sensors. The spatial configuration allowed the measurement of the hydromechanical response of the rock mass at different distances from the shat wall and along different orientations. A velocity survey was carried out between 3 parallel boreholes (REP2301 to REP2303) within the REP experiment, to record initialization and evolution of damage around the shaft during and after the excavation phase.

After the instrumentation of the drift, the shaft sinking resumed. The shaft excavation, with a diameter of 6 m has progressed by blasting. The monitoring focused on 10 different blastholes. They are numbered and presented in Fig. 3. The thickness of the removed rock at each blasthole extends from 1.82 to 3.41 m (28 m 3 to 96 m 3 ) respectively for blasthole 1 and 6. Blasthole 10 (2.76+2.95 m) corresponds to two blastholes, which were reduced to one for the purposes of this study because the acquisition system stopped during this period. For blastholes 2, 6 and 9, the front level was dug deeper with a pneumatic hammer respectively for 40, 70 and 30 cm.

DESIGN AND EMPLACEMENT OF THE EXPERIMENTAL DEVICE

Two key issues have been identified during the design of the velocity device; the way to seal probes with material presenting mechanical properties similar to the claystone and the location of the sensors (distances from the shaft wall, distance between sensors, orientation,…) Some laboratory tests were performed to choose the best ratio of cement and kaolin mixture following a methodology of material reconstitution and ultrasonic measurement under stress [START_REF] Darot | Acoustic wave velocity and permeability, evolution during pressure cycles on a thermally cracked granite[END_REF]. The cement was chosen for its weak ultrasonic damping. The kaolin was chosen for its low swelling effect (for a clay) and especially its capacity to reduce the compressive strength and the mechanical modulus of the mixture. The objective was to attain the mechanical characteristics of the argillite (Young modulus, Poisson ratio) to avoid a discrepancy between the sealing column and the rock media. The ratio of 15% of kaolin was selected as the best compromise between the modulus and the attenuation.

A constitutive model has been proposed by Andra on the basis of Hoek and Brown criterion [START_REF] Souley | Constitutive models for rock mass: numerical implementation, verification and validation in FLAC and Numerical Modeling in Geomechanics[END_REF]. In this generalized Hoek-Brown model, damage in the pre-peak phase is approached by the theory of plasticity and the transition between the brittle failure and the ductile behavior is also integrated as observed from the laboratory characterization on samples. The elastoplastic model uses an associated flow rule with hardening and softening, which exhibits four phases: I-elastic, II-damaged, III-broken, IV-residual. Numerical modelling of the shaft sinking was performed using this model in 3D taking into account the anisotropic in situ stress field and exhibited stresses changes and damaged zones. . On the other hand laboratory tests were performed on samples cored at a depth of -460 m to determine an empirical relationship between mechanical parameters (deviatoric and confining stresses, damaged phase) and the P-wave velocity under different stress path [9][14][16]. This rough empirical law enables the possibility of back calculating, the expected velocity changes around the shaft from the 3D computed stresses in order to determine the location of the probes and sensors.

The experiment device comprises 3 ultrasonic probes 22.41, 22.71 and 23.97 m long (named G1, G2 and G3), installed in 3 parallel and inclined boreholes from the niche at -447 m depth. The target zone of investigation is located near the shaft between the levels -465 m and -475 m. The boreholes were positioned to place the sensors in the zone where the expected velocity changes induced by damage of the rock are the most important. The sensors are located in the last 5 m at the bottom of the probes. The distance between two sensors of the same probe is based on a geometric sequence fitted on the expected velocity change.

The 3 probes comprise 21 ultrasonic sensors and 6 clinometric sensors. The boreholes were drilled end of December 2004, and borehole trajectories were recorded. The velocity survey probes were emplaced using the clinometric sensors to choose the orientation of the ultrasonic sensors along the medians of the device (Fig. 3) . At the end the ultrasonic probes were sealed with the defined mixture of cement and kaolin. The real distance between the ultrasonic sensors and the shaft wall was calculated by linear interpolation between the measured borehole trajectory and the topographic measurement of the shaft wall performed during excavation. This distance ranges from 0.30 to 2.57 m with an accuracy close to 1 cm.

EQUIPMENT

The ultrasonic measurement equipment used in this experiment was developed in the framework of the European OMNIBUS project [START_REF] Pettitt | Development of the Tools and Interpretation Techniques for Ultrasonic Surveys to Monitor the Rock Barrier around Radioactive Waste Packages in Geological Repositories[END_REF]. One unique feature of this equipment is that it is able to switch automatically the sensors from pulser to receiver using customdesigned pulser-amplifier units. This is in spite of the difference in voltages required for the two modes of operation. During the pulsing it is important to deliver the maximum energy into the rock, and in receiving mode high amplification is required as signals are very low because of the strong damping of the rock media. Therefore, the sensors and associated electronics must be adapted to an alternative voltage of 5 V in receiving and 500 V in pulsing.

The OMNIBUS amplifier boosts the ultrasonic signal over a frequency band from 10 to 1000 kHz corresponding to the same frequency band as that of the sensor with a gain from 40 to 60 dB (x100 to x1000). The sensor has been designed to deliver a signal with a frequency peak around 65 kHz corresponding to the optimized frequency for argillaceous media. This frequency is sufficiently high to detect the expected perturbations, but not too high for the energy of the signal to be too damped in a field experiment. The sensors and downhole electronics (pulser-amplifiers) are also designed to withstand the borehole environment, especially regarding resistance to the pore pressures that apply at this depth. They are attached to probes suited to borehole applications. A system with springs provides the coupling of the sensors to the borehole wall while the hooding of the sensors has the same shape as the wall of the borehole.The probes are equipped with systems for injecting sealing material and inclinometric sensors to orientate the body of the probe so that the sensors are in a facing position, so as to maximize the transmitted energy. An interface unit controls the amplifier-sensor pairing. Its first main function is to supply power for the amplifiers and the sensors and to transmit the analogue signals to the acquisition system. Its second main function is to command the running mode of the amplifier-sensor pairing: transmitter or receiver. The other important element of the acquisition unit is the acquisition box that digitizes and transfers the analogue signals to the acquisition computer. 8 GAGE acquisition cards of 2 channels each compose the unit. Each channel samples with a frequency of 10 MHz and a vertical resolution of 16 bits. A computer drives the entire acquisition unit via a software, which controls the continuous monitoring with an automatic acquisition by triggering once a signal threshold is exceeded.

PROCESSING

The P-wave velocity is the main parameter of this experiment; its absolute value is calculated on the basis of the shortest distance between the transmitters and the receivers and the transit time of the first arrival of the full wave. This velocity does not take account of other bent rays created by velocity gradients in the medium. Consequently, the velocity is measured along a straight ray at the peak frequency of the sensor filtered by the rock ( 65 kHz) and is expressed in metres per second. On the same ray path, unlike the absolute velocity, the velocity change is not much dependant on the exact ray path of the wave. Its calculation is based on the cross correlation of the signals for the same pair of sensors as a function of time. The cross correlation is performed on the first phase of the signal and has an uncertainty of around 1 m/s. However, during shaft excavation, the deformation measurements (measured with extensometers) showed deformations of around 1000 to 2000 m/m. The distance between two ultrasonic sensors could therefore change by as much during the experiment. This uncertainty is negligible between two closes in time measurements. However, it represents up to a maximum of 6 m/s from the beginning to the end of the experiment.

The absolute uncertainty of the velocity relative to the ultrasonic layout mainly depends on the uncertainty of the positioning of the sensors, which was estimated at ±5 mm. This uncertainty is maximum for the shortest rays (0.7 m). It was estimated at ±17 m/s in this case. This uncertainty applies when the velocities of different rays have to be compared, but not when one compares velocity variations.

The estimation of the transversal anisotropy of the medium in terms of velocity is based on the hypothesis of a punctual approximation of the device. In other words, velocity variations from one ray to the next only depend on the anisotropy, which is considered as homogeneous over the entire volume. The anisotropy probably changes from one point to another, especially according to the lithology and the deconfinement/damage pairing. However, the number of rays and their spatial distribution are too low, to provide sufficient spatial resolution, so just one anisotropic value is calculated for the entire device. The anisotropy estimation method is based on the non-linear fit of the experimental data following the form function defined by Thomsen [START_REF] Thomsen | Weak elastic anisotropy[END_REF]:

) ( sin ) ( sin ) ( cos ) ( 4 2 2 m m m m B A V V (1)
where:

V P-wave velocity (m/s);

V m maximum velocity (m/s); θ incidence angle of the ray compared to the horizontal (rad); θ m phase angle difference where the velocity is maximum;

A and B coefficients respectively represent the slope and the amplitude of the anisotropic model.

Since the velocity is, in principle, little dependent on the azimuth of the rays, it is possible to assume transversal anisotropy. The anisotropic calculation can therefore be extended to the third dimension on the basis of a symmetrical revolution around the direction of the main minimum velocity, which, in this case, is close to the vertical direction. We sought to generalise the incidence angle of the rays over the whole horizon by assuming:

( θ -θ m ) 2 = ( x -x m ) 2 +( y -y m ) 2 (2)
where x and y are the coordinates of the ray in the stereographic presentation respectively in the east and north direction.

with (1) and (2):

2 2 4 2 2 2 2 2 2 2 ) ( ) ( cos ) ( ) ( sin ) ( ) ( cos ) , ( m m m m m m m y y x x B y y x x y y x x A V y x V (3)
The estimation of the transversal anisotropy is calculated with a Levenberg-Marquart type resolution [START_REF] Levenberg | A Method for the Solution of Certain Problems in Least Squares[END_REF] [START_REF] Marquardt | An Algorithm for Least-Squares Estimation of Nonlinear Parameters[END_REF]. The transversal isotropy model is applied to each series of velocity field measurements, which makes it possible to determine several indicators that take all the data into account. These indicators are determined from the parameters of equation [START_REF] Thury | The Mont Terri rock laboratory, a new international research project in a Mesozoic shale formation[END_REF]. They are defined by the extremes and the phase differences of the equation representing respectively the maximum and minimum amplitude, the azimuth and the dip of the main velocities. The anisotropy model integrates all the ultrasonic rays and extracts the strong trends corresponding to the perturbation affecting the sounded rock volume. The results can be represented graphically by the regression surface that minimizes the discrepancy with the experimental data.

RESULTS

The results presented in this paper correspond to the raw data of 461 measurement runs with 15 to 18 emissions each. That represents about 17,000 different signals over the whole acquisition period from the 13 th January to the 25th July 2005. The total number of transmitted ultrasonic rays was theoretically 272 per run with the final configuration. However, taking into account the distribution of the device over 5 m and the damping of the rock, the total number of rays that could be recorded and interpreted at the start of the acquisition was 115. The maximum distance covered by the ultrasonic wave in the argilites was 1.95 m, consistent with the 2 m maximum observed in the EZ-A and ED-B experiment of the Mont Terri URL [START_REF] Forney | Caractérisation par méthodes ultrasoniques de la zone endommagée induite par le creusement d'un tunnel en milieu argileux : Cas d'étude au tunnel du Mont Terri[END_REF] [START_REF] Damaj | Auscultation et surveillance des perturbations hydromécaniques d'ouvrages souterrains par méthodes ultrasonores[END_REF].

The signal to noise ratio of the recorded signals depends on the orientation and the length of the rays. The noise level was around 70 mV high and the maximal amplitude of the strongest signals experienced distortion and reached 1 V. So the signal to noise ratio varied from 0 dB to 22 dB respectively for 1.95 and 0.76 m ray lengths. The signals were processed in this study for a signal to noise ratio over 5.

x

East North y …θ

NATURAL VELOCITY CHANGES

In the ultrasonic device, some rays constitute 3 by 3 triangles located along the same sub horizontal planes at different depths. The rays between G1 and G2 are very close to parallel to the shaft wall (plane A), the rays between G2 and G3 are close to the perpendicular (plane A') and the rays between G1 and G3 are close to the oblique line ( 45°, plane A''). The first measurements of the P-wave velocities on these rays showed an important change of the original velocity according to depth (fig; 4). For example, the velocity varied from around 300 m/s from -465.15 to -466.19 m (1.05 m) for the sensors G2.2-G3.4 and G2.4-G3.6 respectively. This vertical variability is in accordance with the natural gamma radioactivity, but does not satisfy the hypothesis of homogeneity used for the estimation of the initial anisotropy. This hypothesis of a homogeneous field may introduce a bias in the initial anisotropy estimation. However, this bias is avoided for the analyse of the anisotropy change during the shaft excavation using a differential analysis.

The fitting of the velocity data, according to the ray angle incidence and the previously defined anisotropic model, showed an initial transversal isotropy up to 14.2 % with a phase difference of 17.4° to the East and 3.1° to the North. Therefore, the isotropy is not exactly transversal but presents a slight orientation in the horizontal plane (lithology plane) mainly to the east. This deviation may come from the initial stress field anisotropy. This initial orientation and its consequences are always taken into account in the modelling of anisotropy.

TEMPORAL VELOCITY CHANGES

Fig. 5 shows the global velocity changes on few representative rays between two ultrasonic probes compared to the shaft front level. The P-wave velocity is approximately constant before the excavation of the shaft reaches a level 8 m above the ultrasonic network. A significant velocity decrease occurs as the shaft advance reaches the level of the instrumented network, the magnitude of the velocity decrease measured on the ray paths depending on sensors location. Velocity continues decreasing, moderately and gradually, as the excavation of the shaft progresses deeper. Since the ultrasonic sensors closest to the shaft wall are located between -469.7 and -468.7 m deep, it is therefore logical that the blasthole 6 (between -464.9 and -467.6 m on the 17/05/05 further extended to -468.3 m on the 20/05/05) created the biggest disturbance.

The example of Fig. 6 shows the evolution on 4 typical rays of the impact of the 6 th blasthole. This blasthole caused an acquisition cut where the velocity data are missing. The impact of this blasthole can be broken down into two stages: (1) from A to B, a drastic and generalized decrease of the velocity just after the blasthole and (2) from B to C, an increase for the majority of the rays; only the rays the farthest from the shaft wall continued to decrease or stabilize.

The velocity of the ray G3.8-G2.5 decreases slightly after the 6 th blasthole and stabilizes from the 21/05/05. The velocity of the other rays decreases first, reaches a minimum before starting increasing. The ray G3.4-G2.3 reacts quicker than the others, reaching its minimum velocity value before the 20/05/05 at 4pm, while the two others rays G3.6-G2.3 and G3.4-G2.2 reach it only on the 21/05/05 at 12am. Consequently, the increase in the velocity does not start at the same time for each ray. Fig. 7 presents the stereonets (i.e. velocity changes as a function of ray orientation) of the velocity for different stages of the shaft excavation. The first variations affect the rays with a dip very close to the vertical (shaft front at a depth of -455.9 m and -461.7 m), then at -464.9 m the velocity changes along an axis close to that of the minor horizontal stress (N65°E). This means the rays close to the plane passing through the vertical and with the azimuth N75°E are the most affected. After the most disruptive blasthole between 464.9 and -467.6 m, all the rays are affected by a velocity decrease with a bipolar tendency.

The least disrupted rays were oriented along an axis with an azimuth between N90°E and N130°E and between N270°E et N300°E with, in this same zone, a minimum disruption of rays with the lowest dip;

The most disrupted rays were oriented along the complementary azimuths N60°E and N240°E corresponding to the direction of the minor horizontal stress after the boring of the shaft. These two angular zones of major change tended to extend over subsequent blastholes. The affected rays were the ones with a dip less than 50° and azimuths between N310°E and N80°E and between N130°E and N260°E. A focused analysis of the velocity changes just after the 6th blasthole (Fig. 8) shows that the disturbance is mainly for subhorizontal rays in the direction of the horizontal minor stress, what confirms the general trend of the experiment. The rays in the direction of the major horizontal stress show only a very slight decrease in velocity. The results also show an alignment of the rays most affected by a velocity increase along the N70°E-N250°E direction. At the same time, the rays, in the direction perpendicular to this axis show a slight velocity decrease.

SPATIAL VELOCITY CHANGES

For each ray, the velocity changes are presented according to their distance from the shaft wall (Fig. 9). The distance corresponds to the barycentre of the transmitter-receiver pairing. This spatial presentation has been broken down into three significant steps of the experiment: (1) the initial state corresponds to the average of the velocity for each ray obtained with ten measurements just after a weak velocity increase during the first days of the experiment (shaft front at a depth of -451 m); (2) the intermediate state corresponds to an average obtained for each ray over ten measurement points immediately after the blasthole 6 (shaft front at a depth of -468 m); and (3) the final state based on an average, for each ray, of the last ten measurements of the acquisition (front at a depth of -480 m).

The velocity changes during the experiment were overall largest close to the shaft wall and decreased with increasing distance. Two main zones could be distinguished. A first zone close to the shaft wall, with a majority of rays that showed a velocity change over 30 m/s. This zone extends to a distance of 1.6 m from the shaft wall. The maximum variation was attained at 0.7 m from the shaft wall with -67 m/s. There is a second zone where the velocity varied less than 20 m/s over 1.6 m. The increase in velocity measured on some rays just after the 6th blasthole seems to affect only the sensors in the first metre of rock (zone 1) although the blasthole 6 disturbed rock through 1.6 m (zones 1 and 2 in Fig. 10).

ANISOTROPIC MODEL

The extended transversal anisotropic model explained previously enabled an integration of all the raw data. The model provided the overall evolution of the velocity throughout the analysed rock volume. The results are presented on the curves of the variations of the model parameters as a function of time and shaft front progress.

The curves of the main velocities (minimum and maximum velocities, Fig. 11 and Fig. 12) clearly have a stepped shape with (1) a stable period with initial velocities close to the ones of the unperturbed rock, (2) a drastic decrease during the shaft crossing for the maximum velocity direction and more moderate for the minimum velocity direction and (3) a period with a slight velocity decrease after the shaft has gone past.

The decrease in the velocity is significant but remains moderate, around 40 m/s in the direction of the maximum velocity (lithology plane) and around 30 m/s along the vertical. Therefore, the perturbation of the shaft excavation seems to be more limited along the vertical axis. In addition, the velocity decrease in the lithology plane is mainly concentrated during the period just after the blasthole 6, while the decrease along the vertical is distributed over all the blastholes. The minimum velocity is close to the vertical and is influenced earlier as the shaft front and the associated unloading comes, whereas the maximum velocity, close to the horizontal, is influenced only as the shaft front is going past.

The velocity increase at the beginning of the experiment shows two different flow rates that can be seen on the maximum velocity curve (Fig. 16-a). The first flow rate corresponds to the first days after the installation of the probes with a moderate but rapid increase in the velocity of several metres per second. The amplitude of the second flow rate is similar to the first but the velocity change is slow and spread out over the 100 first days of the experiment.

The evolution of the anisotropic factor (ratio between the minimum and maximum velocities) shows the combined effect of the disturbance of the strongest velocities and the weakest velocities. Fig. 16-c shows a greater sensitivity than the ones provided by the changes in the two main velocities analysed separately, particularly for the minimum velocity. From the first blastholes, the anisotropic factor tends to increase by 0.1%. Nevertheless, the biggest changes occur between the blastholes 3 and 9. During the blastholes 3, 4 and 5, the anisotropic factor increases by 0.3 %, decreases by 1 % during the blasthole 6 and finally increases during the blasthole 7. After this last blasthole, the anisotropic factor decreases to a plateau value of 0.3 % below the initial anisotropic factor.

The main anisotropic directions show some weak but significant changes that drive a rotation to the north. This rotation is more significant in the main direction of the high velocities. The oscillation already observed on the anisotropic factor is also highlighted for its orientation, especially with the dip for the direction of the highest velocity. More precisely, the dip increases while the azimuth decreases before the arrival of the shaft front at the level of the device. When the shaft front goes past the device, then the tendency is reversed to levels higher than the initial state. The relationship between velocity change and shaft position forecasted by numerical modelling 1 m from the shaft wall is similar to the average of the velocity change of the anisotropic model (average of the velocity in all directions). This is highlighted by an initial decrease 10 m behind the shaft front and a typical step after the shaft has gone past. However, the amplitudes of the variations diverge by a factor of 10 (Fig. 12) and are discussed further below.

The small amplitude of velocity variations measured with the velocity survey device is consistent with the measurements performed in the shaft just after the excavation. Indeed, as the shaft front was at a depth of 468.3 m, the excavation was stopped more than 2 weeks to allow different kinds of measurements to be done. Seismic refraction performed on the shaft wall and seismic profiles carried out in a radial borehole showed that velocity variations did not extend deeper than a few tenth of centimetres; the resolution of these methods was lower (the frequency was about ten times lower), but it confirms that velocity variations are very small in the volume investigated by the velocity survey.

DISCUSSION

NATURAL VARIABILITY OF THE VELOCITY

The natural variability of the P-wave velocity in argilites is 12% depending on the depth and 14% depending on the dip of the rays. The natural rock is therefore both heterogenic and anisotropic. This variability is greatly superior to the one induced by the shaft excavation. It is only 2% at the maximum.

The velocity variability as a function of depth also affects the estimation of the velocity variability as a function of the orientation of the rays. The extension of the network needs to be sufficient to have a good angular cover of the rays and so this extension necessarily covers a rock volume characterized by a certain heterogeneity. However, in this experiment, the spatial distribution of the rays was sufficiently random for the angular variability to clearly stand out despite the fact that the angular variability is of the same order of magnitude as the spatial variability. It is this natural variability that is expressed in the form of a mainly transversal anisotropy. It comes from the sedimentary structure of argilites with the series of the lithologic interfaces and the preferential orientation of the mineral layers [START_REF]Référentiel géologique du site de Meuse/Haute-Marne, Tome 4 : Le callovo-oxfordien[END_REF]. The waves propagate quicker along the interfaces and along the mineral layers than when propagating perpendicular.

The average anisotropy estimated by the ultrasonic device between -463 and -469.7 m deep is 14.2%. This is higher than the anisotropy estimated on samples by laboratory tests (10%) at a depth of -460 m in the same facies. This difference can be linked: [START_REF] Martin | Observations of brittle failure around a circular test tunnel[END_REF] to the scale difference between the sample and the in situ rock volume, which is caused by a network of lithological interfaces of size greater than the test sample (70x140 mm), (2) homogeneous damage following the deconfinement and desaturation during the sampling, which decreases the anisotropy by the same amount.

DYNAMIC OF THE VELOCITY CHANGES

Therefore, only 2% of the velocity changes are linked to the deconfinement and the damage to the argilites. The drilling of the boreholes itself probably created a disturbance of the velocity field. This disturbance tends to reach equilibrium through reconfinement within the first few days of the experiment (corresponding to the observed measurement increase) due to the slightly expansive properties of the cement-kaolin mixture and a slow resaturation in the immediate borehole vicinity. The disturbances only caused by shaft excavation (Fig. 13) are mainly characterised by: (1) the initiation of the velocity decrease from the blastholes situated 10 m from the instrumental device; it is greater and earlier on the subvertical rays (in the shaft axis); (2) a maximum of disturbance during the shaft crossing with the higher amplitudes on the rays perpendicular to the shaft wall and for the rays located in the zone between 0 to 1.6 m from the wall; (3) a velocity increase just after the drastic decrease (velocity rebound) around the shaft front up to 1.6 m behind the front and up to 1 m from the wall, especially for the rays oriented in the direction of the minor horizontal stress; and (4) a weak oscillation of the amplitude (1%) and the dip (0.4°) of the anisotropy during the shaft crossing and finally a rotation of the velocity field to the north. These measurements seem to indicate initially a deconfinement of the rock ahead of the front and logically oriented along the vertical stress. It is subsequently followed by a stronger deconfinement and perhaps slight damage with microcracking parallel to the shaft wall. This damage, probably complex due to the heterogeneities of the rock, the irregularities of the shaft wall surface and the corresponding stress field, is nevertheless progressive and inversely proportional to the distance from the shaft wall. It seems to be contained within the first 1.6 m of the rock. Other evidences confirm that damage zone is small, with a low level decreasing with depth. For example, permeability measurements performed in the vincinity of the shaft [START_REF] Su | Hydromechanical coupling phenomena observed during a shaft sinking experiment in a deep argillaceous rock[END_REF] show that the most significant permeability increases (one order of magnitude) are detected in the zone very close to shaft wall up to 1 m, small increases (less than one order of magnitude) remains measurable up to 2 m and to 2.5 m depending the orientation to the stress field. Displacement measurement indicates that displacement in the direction of the major horizontal stress is nearly elastic and extent of plastic zone remains small [START_REF] Su | Hydromechanical coupling phenomena observed during a shaft sinking experiment in a deep argillaceous rock[END_REF]. Numerical analysis (cf. chapter 2.2 and [START_REF] Souley | Modelling of the hydromechanical response of a shaft sinking in a deep claystone[END_REF]) exhibits 1.5 m of damaged zone in the direction of minor horizontal stress. All those data are coherent with the velocity field described herein. The oscillation of the anisotropic factor during the mine-by test can be linked in part to pure elastic phenomena such as reorganization of the stress field and in part to irreversible phenomena such as rock damage in the vicinity of the shaft (induced anisotropy). Laboratory tests on Callovo-Oxfordian argilites [START_REF] Pettitt | Development of the Tools and Interpretation Techniques for Ultrasonic Surveys to Monitor the Rock Barrier around Radioactive Waste Packages in Geological Repositories[END_REF] have shown that the initialization of an anisotropic decrease corresponds to the damage initialization. In fact, the anisotropy in this rock is induced mainly by the asymmetric nature of its layered structure. Inducing a random microcracking tends to reduce the difference between the velocity along the layers and across the same layers. This phenomenon is further intensified by the preferential orientation of the deconfinement perpendicular to the shaft wall and the lithology.

The reorientation of the velocity field is due to the rotation of the stress field during the shaft crossing. The minor horizontal stress h tends to rotate towards the shaft axis, in which the radial stress on the wall surface is zero (assuming the argilites are not very plastic). At the same time, the major horizontal stress H tends to rotate towards an axis orthoradial to the shaft with a maximum close to the surface. The velocity changes depend on the value of the ratio of the applied stress divided by the tensile strength of the rock [START_REF] Lockner | Changes in seismic velocity and attenuation during deformation of granite[END_REF] and the direction of propagation of the studied waves. In fact, depending on the level of the applied stress, either the existing cracks will close (low stress), or new cracks will be created (high stress). In the case of weak stress, it is the cracks perpendicular to said stress that will close preferentially. In the case of high stress, cracks parallel to said stress are created, which induces a decrease in the velocity of the waves that propagate perpendicularly to the stress direction [START_REF] Nur | Stress-induced velocity anisotropy in rocks: an experimental study[END_REF].

The rock is clearly deconfined along the radial axis of the shaft and the deviatoric stress could create some microcracking parallel to the shaft wall. However, it is very difficult to determine the contribution of each phenomenon on the wave velocity decrease. The laboratory tests and the numerical modelling used to forecast the velocity changes have overestimated their amplitude by an order of magnitude. This difference can stem from the stress field modelling, especially as regards the uncertainty of the input parameters used at that time (modulus level, stress anisotropy, etc.). The difference also can stem from the laboratory tests, which are not sufficiently relevant to represent damage phenomena on a large scale. Triaxial tests that exactly reconstitute the same stress path in the sounded rock volume would doubtless be more representative. In addition, the duration of the tests could be the determining factor of this difference. The laboratory tests were performed over several hours, whereas the mine-by test lasted several weeks. Different phenomena such as the viscous behaviour of argilites could reduce the initialization or the opening of the cracks. Whatever the case, the velocity changes are slight and consequently the damage generated by the shaft is low. Furthermore, the instrumental device is at the same time too big to sound a layer of homogeneous argilites and too small to sound the whole rock medium around the shaft. The results of this experiment focused on a limited area next to the shaft, so it is not yet possible to know if the extent and the level of the perturbation vary around the entire shaft and, in particular, if the perturbation is lower on the side of the direction of the major horizontal stress.

VELOCITY REBOUND

A non-negligible transitory phenomenon seems to be associated with the irreversible damage and stress release phenomena in the immediate vicinity of the front (1 m). This is reflected by an increase in velocity and seems to propagate with, for example, an increase of 7 m/s for example for the ray G3.4-G2.2. This phenomenon is clearly highlighted by the change of the velocity anisotropy. However, this phenomenon is difficult to understand at the state of our present knowledge and simply on the basis of the ultrasonic velocity. However, we can put forward three hypotheses. This increase could be due to: (1) changes in pore pressure; in fact in the direction of the horizontal minor stress instaneous over pressure are observed due to hydromechanical coupling under an anisotropic state of stress with a low permeable media [START_REF] Su | Hydromechanical coupling phenomena observed during a shaft sinking experiment in a deep argillaceous rock[END_REF]. Over pressure are observed before the front reaches the level of the sensor and when the front crosses the level of a sensor, high peak of pressure appears followed by a pressure drop down. [START_REF] Pettitt | The assessment of Damage Around Critical Engineering Structures Using Induced Seismicity and Ultrasonic Techniques[END_REF] The displacement of the stress front ahead of the shaft front during the progress of the shaft by punctual attack with pneumatic hammer; however the punctual attack with a pneumatic hammer only concerns the last meter of the velocity survey device, while the rebound is observed on rays at the level or even above the shaft front. (3) Blast effect which induces stress redistribution around the measurement boreholes which can affect the measurements. After, the creep behaviour could explain the increase of velocity.

As the next excavation steps resume it is difficult to say if the rebond is just a way to reach equilibrium or if it is linked with a transitory phenomemum. Further information and analysis are required to better explain it.

CONCLUSIONS

The excavation of the main shaft in the argilites of the Meuse/Haute-Marne URL significantly disrupted the velocity field. Monitoring of the P-wave anisotropy during the mine-by test enabled the singular behaviour of the ultrasonic rays to be eliminated while retaining information from each one.

The velocity field deteriorated from 10 m ahead of the shaft front with a maximum disruption during the shaft crossing. This disturbance decreased progressively up to 1.6 m from the shaft wall. The shaft crossing creates an oscillation of the main velocity directions with a significant velocity drop initially and then a differed velocity rebound. Finally, the order of magnitude of the velocity changes showed that deconfinement and damage to the rock were limited. It confirms that the level of damage is low and that the EDZ extent remains small. Those results are coherent with other measurements performed in the vicinity of the shaft (permeability and deformation measurements and other seismic methods) and with elastoplastic 3D modelling (even if measured velocity field is not in accordance in term of amplitude with the one back calculated from the 3D modelling results with an empirical relationship).

The experimental device is very sensitive but cannot clearly separate the deconfinement and damage contribution in the disruption and only sounded a limited area around the shaft. The measurement of the damping of the waves, the transmission of a full shear wave signal and the integration of an orthotropic component in the anisotropic model would improve understanding of the microcracking and the rock damage phenomena.

Figure 1 :

 1 Figure 1: Presentation of the ANDRA URL.

TimeFigure 3 :

 3 Figure 3: Depth changes of the shaft front during the mine-by test.

Figure 3 :

 3 Figure 3: Final experimental device with recorded ray paths and orientation of the probes (view perpendicular to the borehole axis).

Figure 4 :

 4 Figure 4: P-wave velocity changes for horizontal rays in comparison with natural gamma radioactivity (data provided by ANDRA and recorded by SEMM).

Figure 5 :

 5 Figure 5: Temporal velocity changes on few representative ray path during the mine-by-test.

Figure 6 :

 6 Figure 6: Effects on velocity during the 6 th blasthole.

Figure 7 :

 7 Figure 7: Stereonets of the velocity changes at different stages of the shaft excavation.

Figure 8 :

 8 Figure 8: Stereonets of the P-wave velocity just after the 6th blasthole on the right (B-A) and for the velocity rebound on the left (C-B).

Figure 9 :

 9 Figure 9: P-wave velocity change according to shaft wall distance for the initial, intermediate and final states.

Figure 10 :

 10 Figure 10: Velocity changes as a function of the distance from the shaft wall during the 6 th blasthole.

Fig. 11 :

 11 Fig. 11: (a) P-wave change of the maximum velocity versus time (on the left) and shaft depth (on the right). (b) P-wave change of the minimum velocity versus time (on the left) and shaft depth (on the right) (c) Anisotropic factor change versus time (on the left) and shaft depth (on the right).

Fig. 12 :

 12 Fig. 12: Comparison between the expected velocity change at 1 m from the shaft wall and the average velocity change observed on the device as a function of the shaft front progress.

Fig. 13 :

 13 Fig. 13: Interpretative figure of the velocity changes during the mine-by test.


Laura Chiappini
email: laura.chiappini@ineris.fr

Laure Delery

Eva Leoz

Bruno Brouard

Yoann Fagault

A first French assessment of population exposure to tetrachloroethylene from small dry cleaning facilities

Keywords: cleaning facilities, indoor air quality, passive sampling, population exposure

Used as a solvent in the dry cleaning industry, tetrachloroethylene can be a pollutant of residential indoor air, which can cause long term harmful exposures because of its neurotoxicity and probable carcinogenicity. In France, dry-cleaning facilities are integrated in urban environments (shopping malls, residential buildings) and can contribute to tetrachloroethylene exposure for customers and residents. This exploratory work presents the results from five studies carried out in one shopping mall and four residential buildings housing a dry cleaning facility. These studies involved dry cleaning machines fitted with a Carbon Adsorber (CA) and unfitted , with or without Air Exhaust System (AES). Samples were collected in the cleaning facilities and in the apartments located above with passive samplers allowing measurement of time-integrated concentrations on a 7 days sampling period. It has obviously shown the degradation of indoor air quality in these environments and underlined the contributing role of the machine technology and ventilation system on the amount of released tetrachloroethylene in the indoor air. To temper these results, it must be pointed out that some parameters (building insulation, amount of solvent used…) which would influence tetrachloroethylene fugitive release have not been quantified and should be looked at in further studies.

Key words: tetrachloroethylene, dry

Introduction

Tetrachloroethylene (C 2 Cl 4 ), also called perchloroethylene (PCE), is a polychlorinated solvent abundantly used in French dry cleaning facilities (Poirot et al., 2007). It has been identified as an hazardous urban air pollutant causing acute, chronic and potential carcinogenic health effects by the U.S. Environmental Protection Agency (USEPA, 2000). Because of its high volatility, inhalation is the main human route of exposure to PCE. An international target indoor air concentration of tetrachloroethylene has not been set yet. The New York State Department of Health (NYSDOH) has developed a health-based guideline of 100 µg m -3 for chronic exposure (15 ppb) (NYSDOH, 1997). In Germany the value of 5 µg m -3 has been proposed for chronic exposure [START_REF] Schleibinger | Recommandations for establishing target values and guidance values for volatile organic compounds (VOC) in indoor air[END_REF]. In France, there were close to 7500 dry cleaning establishments in 2002 operating 1 to 2 machines. Because of the lack of national data, a complete inventory of the French dry cleaning machine number and type cannot be done. However, it can be assumed that only a few percent of these machines are equipped with a carbon adsorber, CA, to capture PCE vapours, and also a few percent of the dry cleaners are equipped with an Air Exhaust System, AES, (Déléry, 2007). Moreover, in France, dry-cleaning facilities are very frequently located in environments such as large shopping malls or residential buildings. Thus, it seems reasonable to suppose that dry cleaning facilities can significantly contribute to PCE exposure for supermarket customers and building residents. Indeed, when the mean indoor PCE concentration measured in France is 7.4 µg m -3 (Observatory on Indoor Air Quality (OQAI), 2006), PCE levels in residential air can be much higher and range from a few hundred to thousands of microgram per cubic meters [START_REF] Altmann | Neurobihavioural and neurophysiological outcome of chronic low-level tretrachloroethylene exposure measured in neighborhoods of dry cleaning shops[END_REF][START_REF] Delery | Etude du parc français des installations de nettoyage à sec suite à la mise en oeuvre de l'arrêtétype 2345 (2 mai 2002)[END_REF][START_REF] Schreiber | Apartment Residents and day care workers exposures to tetrachloroethylene and deficits in visual contrast sensitivity[END_REF]. Because dry-cleaning can significantly contribute to tetrachloroethylene exposure of neighbourhoods, the aim of this work was to determine indoor air levels in typical French environments housing a dry cleaning facility. Measurements were carried out in one shopping mall and 4 residential buildings with on-site dry cleaners located in the basement floor. The influence of the dry cleaners technology (equipped or unequipped with a CA) and shop ventilation was also assessed during this study.

Study design

Indoor air PCE levels from four residential buildings and one shopping mall with on-site dry cleaners were determined between 2002 and 2007. The sites were chosen to study the exposition of population to tetrachloroethylene in the two main French environment in which dry cleaning facilities are settled: residential buildings and shopping mall. Therefore, two preliminary sampling campaigns using canisters (10 hours sampling period) were conducted in a supermarket (study 1) and in a residential building (study 2) in 2002. An other criterion to chose the sampling sites was the machine technology: in each case, the machine were closed with refrigeration cooling and were equipped or unequipped with a Carbon Adsorber (CA). Thus, for the remaining studies (studies 3, 4 and 5), led between 2005 and 2007, sampling campaigns were performed in three different residential buildings with on-site dry cleaners equipped (study 5) and unequipped (studies 3 and 4) with a CA. During these campaigns, 7 days time integrated samples were collected using passive samplers in order to obtain more reliable estimation of residents exposure.

The sites description as well as four important parameters (machine technology, dry cleaning facility ventilation, shop activity and temperatures in the shop and the apartments, when available) influencing the amount of fugitive PCE concentrations in indoor air are reported in Table 1. In each case, only one dry cleaner was in function in the facility. Studies 1 and 5 concerned dry cleaners using CA equipped machines. Air exhaust system was only present in study 5. In study 1, samples were taken in 10 different locations of a supermarket in a shopping centre housing a dry cleaning facility. Studies 2, 3 and 4 concerned dry cleaners using unequipped CA machines. The dry cleaning activity was high in study 1 (11 loads/day), normal in studies 3, 4 and 5 (6-8 loads/day), low in study 2 (3 loads/day). No measurement was performed in the dry cleaning shop of study 1.

Sampling was only performed in the staircase of the residential building in study 2 while sampling in study 3, 4 and 5 was performed in both staircase and apartments.

No background site was chosen to give information on background PCE levels since it has been assumed that the dry cleaner was the major source of PCE. In fact, it must be pointed out that, for each chosen dry cleaning facility, no other dry cleaner, closed enough to affect PCE measured level, was present in the surrounding area. Moreover, all of them were located in residential area, far from potential industrial influence. Mean PCE concentrations measured in the dwellings of 567 French families between 2003 and 2005 (Observatory on Indoor Air Quality (OQAI), 2006) will be taken as a reference to compare the studies.

Air monitoring and sample analysis

Sampling devices and period for each study are shown in Table 2. Studies 1 and 2 were carried out with 6 L stainless steal canisters while studies 3, 4 and 5 were carried out with passive samplers. For each studies, all the samples were collected at the same time.

Sampling equipment and conditioning

Electropolished stainless steel canisters were purchased from Resteck (Silocan TM ). Cleaning was performed by repetitive evacuation and refilling with humidified zero air at 100 °C. The efficiency of the cleaning procedure was verified by GC analysis. A full description of the adopted procedure can be found in [START_REF] Castellnou | Refrigerated multibed adsorption in sampling and analysis of atmospheric light hydrocarbons at ppb (v/v) and sub-ppb (v/v) concentrations[END_REF] Passive tubes were conditioned by purging with ultra-high-purity helium at 50 mL min -1 while heated at 350 °C for at least 8 hours. After conditioning, the sorbent tubes were capped with stainless steel Swagelock nuts.

Sampling and analysis with canisters

Canisters were fitted with Veriflow controller (Restek inc) maintaining a constant flow (about 10 mL min -1 ) into the canister over the sampling period of 10 hours. Canisters were analysed by transferring an aliquot (2 litres) of the collected air onto a Perkin Elmer desorption tube filled with 150 mg of carbotrap B and cooled down to -10°C. Details of the transfer apparatus were described by [START_REF] Castellnou | Refrigerated multibed adsorption in sampling and analysis of atmospheric light hydrocarbons at ppb (v/v) and sub-ppb (v/v) concentrations[END_REF] Samples tubes were then subsequently analysed with a Perkin Elmer automated thermal desorber ATD 400 coupled with a gas chromatograph equipped with a flame ionisation detector (FID). Analytical conditions are given in Table 3.

For this sampling time (10 hours) and flow rate (10 mL min -1 ), the quantification limit of this method is 4 µg m -3 of tetrachloroethylene. The uncertainty associated to this measurement is estimated to 30 %.

Sampling and analysis with passive tubes

Silent and not cumbersome, this very convenient technique for indoor air measurements allows time-integrated concentrations over long sampling periods. A key factor when performing measurements with passive samplers is the steadiness of the uptake rate over the sampling period. The performance of a radial sampler consisting of the combination of an empty Perkin Elmer tube fitted with a 4.8 mm od Radiello cartridge containing 150 mg Carbotrap B obtained from the Fondazione Salvatore Maugeri (Radiello Perkin-Elmer sampler) and a radial passive sampler (Radiello model 145) was evaluated in a 150 L Pyrex exposure chamber under controlled conditions of temperature, relative humidity and wind speed. A full description of the exposure chamber is given in [START_REF] Gonzalez-Flesca | A new laboratory test chamber for the determination of diffusive sampler uptake rates[END_REF]. The chamber operating conditions were 20°C, 50 % relative humidity, 1 m s -1 wind speed and the VOCs concentrations were 85 µg m -3 of benzene, 132 µg m -3 of trichloroethylene, 180 µg m -3 of toluene and 127 µg m -3 of tetrachloroethylene. As schematic of the tubes is given on Figure 1.

The uptake rate steadiness for both samplers, Radiello 145 filled with carbograph 4 used with the yellow diffusion membrane (polyethylene 5 mm thick, 10 ± 2 µm porosity, 150 mm diffusive path length), and Radiello-Perkin-Elmer (also filled with carbograph 4) and several compounds (benzene, toluene, trichloroethylene and tetrachloroethylene) as a function of sampling time is presented in Figure 2. It can be observed that the axial Radiello Perkin-Elmer sampler provides a lower but steadier uptake rate than the radial sampler. The sampling flow rate for PCE has then been determined equal to 0.3 mL min -1 . No experience has been carried out in the chamber at high PCE concentrations (about a few mg m -3 ) to estimate the tube limit of linearity range. However, some field experiments in a dry cleaning facility involving active sampling on tubes filled with Carbopack X (CPX) adsorbent and Radiello-Perkin-Elmer passive tubes have been performed. Passive samplers have been exposed during 7 days while several samplings, 9 series of 2 samplings covering the 7 days passive tube exposure, have been performed. Both methods gave the same results (93 mg m-3 has been measured by active sampling while 97 mg m -3 has been measured by passive sampling). It can thus be concluded that Radiello-Perkin-Elmer tubes can be used for sampling in high concentration PCE atmosphere. This sampler was therefore chosen for field measurement campaigns. Samples were thermodesorbed by an ATD 400 Perkin-Elmer thermodesorber, pre-concentrated on a cold trap (Turbomatrix Perkin Elmer pre-concentrator) and analysed by GC-FID (Autosystem XL Perkin Elmer equipped with a capillary column PLOT Alumine /KCl) For this sapling time (7 days) and flow rate (0.3 mL min -1 ), the quantification limit of this method is 0.1 µg m -3 . The uncertainty associated to these passive tube measurements is estimated to 30 % (the uncertainty was determined on the bases of 3 chamber experiments performed in the same conditions as discribed before and on the bases of the comparison between active sampling as reference method and Radiello-Perkin-Elmer passive sampling during the field experiment described above).

Calibration and quantification

Quantification was performed by using PCE liquid solutions. 1 µL of each standard solution were evaporated on Radiello-Perkin-Elmer tubes with a Dynatherm tube conditioner Model 60 (Supelco). PCE vaporised mass were ranging from 3 to 130 µg corresponding to 0.9 to 45 mg m -3 As for measurement of lower PCE concentrations (about a few µg m -3 ), quantification was performed by on-line injection of hydrocarbon gas standards (50 ppb) mixture from Air Liquid (CARBAGAS, Switzerland).

Results and discussion

A summary of the results obtained during the five studies is presented in Table 4. Standard deviation is given for each measurement except for studies 1 and 2 since only one canister sampling has been made at each sampling point. Also, only one measurement has been made on the building floors during study 4. Lets remind that the measurement uncertainty including sampling and analysis is estimated to 30 %.

Tetrachloroethylene concentrations in the shopping centre (Study 1)

The highest tetrachloroethylene concentration is naturally measured in front of the dry cleaning shop (678 µg m 3 ). Relatively high concentrations can still be measured at the sampling points far away ( the furthest sampling point is about 60 m away from the dry cleaning facility) from the emission source (50 µg m 3 at the back of the supermarket). Figure 3 presents an isoconcentration plot calculated with Surfer (Scientific Software Group) in order to estimate the concentration profile in the supermarket. Despite the use of a machine equipped with a CA, the fugitive emission of PCE is not negligible. A suitable AES is recommended to further reduce the emission of PCE and thus the exposure inside the shopping mall. Indoor air concentrations in shopping mall housing dry cleaning facilities is not documented. Concentrations levels measured in this study should encourage further research.

Tetrachloroethylene concentrations in the dry cleaning shops (Studies 2-5)

Several measurements were carried out in each dry cleaning facility (at least 3 samples per shop). The tetrachloroethylene arithmetic mean levels are reported in Table 4.

Except study 2 (low dry cleaning activity, special sampling period and technique), studies 3-5 concentrations ranges indicate rather low levels of PCE in dry cleaning shops (8 to 53 mg m -3 ) compared with occupational literature data. Occupational mean concentrations levels of 69 mg m -3 (10 ppm) were reported by Poirot et al., 2007 in a 3 dry cleaners study.

In studies 3 and 4, PCE mean concentration levels are similar considering measurement uncertainties (32 10 and 53 16 respectively). However, PCE concentrations measured in study 4 are slightly higher than concentrations measured in study 3, which could be explained by the higher temperature in study 4 (26.9 °C in the dry cleaner and 26.9 in the apartments) than in study 3 (19.2 °C in the dry cleaner and 18 in the apartments). Both of them were in the same configuration: no AES, CA unequipped machine, normal activity. However, significant decrease of tetrachloroethylene concentrations in dry cleaning facilities operating with a CA equipped machine (study 5) in comparison to an unequipped machine (study 3 and 4) can be observed even if other contributing parameters (machine life span, machine maintenance and PCE consumption) have not been assessed. These results enlighten the positive influence of a CA and an AES on indoor air quality.

It can be noticed that measured concentrations in study 2 also carried out in a dry cleaner involving an unequipped machine are comparable to concentrations measured in study 5 (4 th generation machine). This could be explained by the dry cleaning facility activity during the sampling period which was low in the case of study 2 and normal in the case of study 5. These two studies are difficult to compare since the sampling periods were different, only 10 hours in the case of study 2 and a whole week in the case of study 5. Samples in study 2 were not as representative of population exposure as samples in study 5 and may have not been taken when PCE concentrations in the dry cleaner were the highest. Besides, it should be noted that dry cleaners PCE levels measured in this study were far below the 335 mg m -3 exposure limit value required by the French occupational health regulation.

Tetrachloroethylene concentrations in the residential buildings (Study 2-5)

The concentration ranges measured in the apartments and on the floors from the bottom up to the top of the buildings are reported in Table 4. A decrease of PCE concentrations with floor levels can be observed as shown on Figure 4 and Figure 5. Mean concentrations are represented on each figure with standard deviation which is about 5 % (ranging from 1.5 to 6.4 %) for floors measurements and 10 % (ranging from 1 to 12 %) for apartment measurements. Mean values, 5 and 10 % are taken for graphic representation. The lower floor standard deviation can be explained by the fact that replicates measurements on the floor are taken exactly at the same location whereas in the apartments some sampling are made in different places.

Comparing concentrations measured inside the apartments and on the building floors, it appears that, in each case (studies 3-5), concentrations are greater on the floors than inside except for level 1 which shows concentrations 20 to 50 % greater inside than on the floors. This could be explained by the fact that level 1, located directly above the dry cleaning facility, is the level under the greatest PCE influence and that accumulation phenomena are more important in apartments than in the staircase.

Significant lower indoor air levels in both apartments and dry cleaners can be observed in study 5 (CA equipped machine) in comparison with studies 3 and 4 (CA unequipped machine). Indeed, mean concentrations measured in the apartment located directly above the dry cleaning facility ranged from 296 µg m -3 in the case of a CA equipped machine to 2.9 mg m -3 in the case of an unequipped machine. This reduction could be attributed to the combined action of the CA and the AES. Residential indoor air levels found in dry cleaner buildings with type unequipped machines are consistent with published value. Indoor air concentrations up to 5.5 mg.m 3 were reported in Germany, New Jersey and New-York respectively by [START_REF] Altmann | Neurobihavioural and neurophysiological outcome of chronic low-level tretrachloroethylene exposure measured in neighborhoods of dry cleaning shops[END_REF][START_REF] Delery | Etude du parc français des installations de nettoyage à sec suite à la mise en oeuvre de l'arrêtétype 2345 (2 mai 2002)[END_REF][START_REF] Schreiber | Apartment Residents and day care workers exposures to tetrachloroethylene and deficits in visual contrast sensitivity[END_REF] In all cases (studies 2 to 5), these values are higher than the mean PCE concentrations 7.4 µg m -3 measured in the dwellings of 567 French families between 2003 and 2005 (Observatory on Indoor Air Quality (OQAI), 2006).

Discussion: measurement uncertainties

The reduction of PCE concentrations noted in study 5 can be attributed to the use of a carbon adsorber and mechanical ventilation. Nevertheless, it must be pointed out that indoor air PCE levels can depend on several other parameters that have not been quantified in this study such as the air renewal rate not only of the dry cleaning facility but also the apartment located above, the residents life habits (smoking habit, ventilation…) which have only been evaluated by a simple questionnaire given to the residents, the building structural conditions (e.g. insulation, poorly sealed pipe chases or crack in the wall) the amount of solvent used,… all these factors should be taken into account for a full understanding of PCE fugitive release process in dry cleaner buildings.

Conclusion

These five exploratory studies have shown the degradation of indoor air quality in one shopping mall and four residential buildings housing a dry cleaning facility. Tetrachloroethylene concentrations in residences in the same building were found to decrease with floor levels. The results have underlined the contributing role of the dry cleaning machine technology and ventilation on the amount of tetrachloroethylene fugitive emissions. Among the 4 dry cleaner building studied, the lowest indoor air levels were measured in the residential building with onsite dry cleaners equipped with a CA fitted machine and a mechanical air exhaust system. For the three other dry cleaner buildings equipped with CA unfitted machines, two of them had indoor air concentrations at all floor levels (5 and 7 floors buildings) all above 1 mg m 3 . Such concentrations levels are much higher than the health-based guideline of 100 µg m -3 (15 ppb) recommended by the New York State Department of Health (NYSDOH). Although this study was performed on a limited number of dry cleaning establishments, the finding confirmed that indoor air concentrations in some residential dry cleaner buildings are of public health concern and that measures should be taken to reduce the exposure. Efficient ventilation and the use of CA fitted machine tend to minimise the exposure of the occupant in building with on-site dry cleaners. However, it should be pointed out that the impact of each factor has not been quantified separately. In addition, other potential important parameters such as building insulation or maintenance practice on the amount of PCE available for diffusion through the building have not been quantified at all. Further investigation is needed to identify factors affecting PCE concentrations in order to take appropriate measures to reduce the exposure of dry cleaner building residents. A population risk exposure based the result from this study will be presented in a separate article. This is a first French study on the population exposure to dry cleaning tetrachloroethylene as a function of machine technology. Supported by the French Ministry of Environment, this study

Figure 1 :Figure 2 :Figure 3 :

 123 Figure 1 : Schematic of the Radiello Perkin Elmer configuration (Gonzalez-Flesca and Frézier, 2005)

Figure 4 :

 4 Figure 4 : Concentrations measured on the floors from bottom to top of the residential buildings (studies 2 to 5). Standard deviation for passive tube measurements (studies 3-5) are given (no standard deviation is given for studies 2 and 4 since 1 sample per location has been taken).

Acknowledgement

The authors would like to thank the French Ministry of the Environment which have supported this work

has important insight into the issue concerning the possibly harmful health effects of the population indoor exposure to tetrachloroethylene.


Frédéric Y Bois
email: frederic.bois@ineris.fr

Systems biology GNU MCSim: Bayesian statistical inference for SBML-coded systems biology models

Statistical inference about the parameter values of complex models, such as the ones routinely developed in systems biology, is efficiently performed through Bayesian numerical techniques. In that framework, prior information and multiple levels of uncertainty can be seamlessly integrated. GNU MCSim was precisely developed to achieve those aims, in a general nonlinear differential context. Starting with version 5.3.0, GNU MCSim reads in and simulates SBML models. Markov chain Monte Carlo simulations can be used to generate samples from the joint posterior distribution of the model parameters, given a dataset and prior distributions. Hierarchical statistical models can be used. Optimal design of experiments can also be investigated.

INTRODUCTION

Most systems biology models focus on simulating the behavior of pathways or networks. Studying the properties of theoretical systems is a worthy endeavor, but uncovering the characteristics of actual systems, or testing hypotheses, requires a step beyond pure simulation, into statistical inference. Such inference is needed to correctly account for most sources of uncertainty that stem from experimentation (dealing with uncertainty about model structure is somewhat more difficult but can be partially tackled through model checking.) Inference about the parameter values of complex models, well represented in systems biology by large sets of differential equations, is efficiently performed with Bayesian numerical techniques [START_REF] Gelman | Physiological pharmacokinetic analysis using population modeling and informative prior distributions[END_REF][START_REF] Rogers | Bayesian model-based inference of transcription factor activity[END_REF], in which prior information and multiple levels of uncertainty can be seamlessly integrated. GNU MCSim is a numerical simulation and Bayesian statistical inference tool for algebraic or differential equation systems [START_REF] Bois | MCSim: A simulation program[END_REF]. Other programs have been created to that end: Matlab is one of the most general and easy to use, and BioBayes (Vyshemirky and Girolami, 2008) is a noteworthy specialized software. Still, there is room for free GPL software able to efficiently run computer intensive multilevel Bayesian analyses. GNU MCSim was created specifically to perform Monte Carlo analyses in an optimized, and easy to maintain environment.

PROGRAM OVERVIEW

GNU MCSim is written in ANSI-standard C language. The source code is freely available under GNU GPL and can be compiled for any system, provided an ANSI C compliant compiler. The code can be freely modified.

The software consists in two parts. A model generator, "mod", takes a user defined model written in SBML or in GNU MCSim native syntax, and translates it into C. A dynamic library of routines is then linked to the C version of the model. The resulting executable can run Monte Carlo simulations under a variety of conditions, using a specified statistical model.

Defining structural models in native syntax

The native syntax allows you to specify a set of linear or nonlinear algebraic equations or ordinary differential equations to solve. The model file in that case is an ASCII text file with several sections, including global declarations of variables and parameters, dynamics' specifications (with derivative calculations), model initialization, and output computations.

2

Reading SBML models and applying a template

The Systems Biology Markup Language (SBML) provides a standard for representing biochemical pathway models [START_REF] Hucka | The systems biology markup language (SBML): a medium for representation and exchange of biochemical network models[END_REF]. It can code for models of metabolism, cell-signaling, transcription etc. SBML is maintained since the mid-2000 by an international community of software developers and users (cf. http://sbml.org).

GNU MCSim reads SBML Level 1 models and most Level 2 models. Some features of Level 2 are not yet recognized by mod, such as "func-tionDefinition", "unit", and "rateRule"). Compartments are ignored unless a model template for circulating species is given.

It is possible to read in several SBML models at once. In that case, they are merged, in the sense that if a chemical species is present in several of them, its rate equation is constructed from the reaction descriptions of each SBML model. Obviously, in that case, the concurrent models should be coherent in their structure and notation.

Performing simulations and statistical inference

After compilation of a model into a GNU MCSim executable, five types of simulations can be run, as specified in simulation input files:

Simple runs produce time-courses of the model variables. The main integration routine is Lsodes, from the SLAC Fortran library [START_REF] Gear | The automatic integration of ordinary differential equations[END_REF].

Monte Carlo simulations perform repeated runs across a randomly sampled region of the model parameter space.

Markov-chain Monte Carlo (MCMC) simulations simulate a Markov chain in the model parameter space [START_REF] Gilks | Markov Chain Monte Carlo In Practice[END_REF]. They can provide samples from a Bayesian joint posterior distribution of the parameters, given a structural model (e.g. an SBML model), a statistical model (specified in a simulation file), prior parameter distributions and data for which a likelihood function can be computed. Fast Metropolis (vector) sampling is available, as well as Metropolis within Gibbs sampling (one component at a time). GNU MCSim handles multilevel (e.g., random effects and mixed effects) statistical models [START_REF] Gelman | Physiological pharmacokinetic analysis using population modeling and informative prior distributions[END_REF], particularly useful when several datasets are at hand, or when variance structures are complex.

Set-points simulations solve the model for a series of specified parameter sets, placed on a regular grid, for example, or generated from previous Monte Carlo or MCMC simulations.

The Optimal design procedure optimizes the number and location of observation times for experimental conditions, minimizing the variance of a set of parameters or outputs, given a structural model, a statistical model, and prior parameter distributions [START_REF] Bois | Optimal design for a study of butadiene toxicokinetics in humans[END_REF][START_REF] Amzal | Bayesian optimal design via interacting MCMC[END_REF].

Application examples

The yeast glycolysis model (YeastGlycolysisJDClean.xml) of Pritchard and Kell (2002), provided with JDesigner 2.1c and the Systems Biology Workbench 2.7.8 (http://sbw.kgi.edu/index.htm), was used. Fig. 1 shows the model (17 state variables and 95 parameters) and the results of time-course simulations over a minute. Those simulation results are identical to those obtained with the Jarnac simulator of SBW 2.7.8. To generate a synthetic dataset, Gaussian error (10% CV) was added to points along each species trajectories at 7 times (including time 0.) Those points, used as "data" in the subsequent analyses, are also plotted on Fig. 1. Posterior distributions were obtained by MCMC sampling for 10 model parameters (see details in the supplementary material.) Independent uniform priors and lognormal likelihood functions were used. Fig. 2 shows results obtained with 3 concurrent chains, in two steps each (10000 runs with component sampling, 1 million runs with vector sampling.) Vector sampling was more efficient, given the very strong posterior covariance between some model parameters, for example between J0_Vmax and J0_kglc. Convergence was achieved, with all R diagnostics [START_REF] Gelman | Physiological pharmacokinetic analysis using population modeling and informative prior distributions[END_REF] equal to 1, over the last 5 10 5 iterations. Each chain took about 80 minutes on a i686 machine clocked at 3.6 GHz. Keeping one sample for every 50 in the last 5 10 5 of the 3 chains led to a joint posterior sample of 30000 values for each parameter, whose empirical percentiles are shown Fig. 2. Due to strong posterior covariances, the posterior samples tend to be higher than the true parameter values used to generate the simulated data. An application of design optimization to fix that problem is given in the supplementary material, together with an example of multilevel modeling.

CONCLUSIONS

GNU MCSim can efficiently perform Bayesian inference and optimal design computations for systems biology models. Multilevel models can also be used. Despite a few limitations in reading in SBML level 2, which will be removed in a coming version, GNU MCSim can import most SBML models and allows users to take advantage of its speed and flexibility. Its code is open and can be freely tailored to specific needs.

Fig. 1 .

 1 Fig. 1. JDesigner representation of the yeast glycolysis model of Pritchard and Kell (2002) (left). On the right is a simulation of the time course of all model species (note the logarithmic x-axis.) The dots mark the values taken as simulated data for the Bayesian calibration exercise.

Fig. 2 .

 2 Fig. 2. Trajectory of 3 Markov chains for parameter J0_Vmax at convergence over the last half-million runs (left); On the right: percentile plot (5 th , 25 th , 50 th , 75 th and 95 th percentiles) of the posterior sample obtained for each parameter, overlaid with their "true" value (crosses).

ACKNOWLEDGEMENTS

Funding: This work was supported by the European Commission 6 th Framework Program, Priority 6 (Global change and ecosystems), project 2-FUN [contract #036976]; and the French Ministry of Sustainable Development Program [BCRD 2004 DRC05].

The GNU GPL source is available at http://savannah.gnu.org/projects/mcsim. A distribution package is at http://www.gnu.org/software/mcsim. GNU MCSim is written in standard C and runs on any platform supporting a C compiler.