J Sarout
email: sarout@geologie.ens.fr.

L Molez

Y Guéguen

N Hoteit

Shale Dynamic Properties and Anisotropy under Triaxial Loading: Experimental and Theoretical Investigations

Keywords: shale, elastic wave velocity, anisotropy, damage, effective medium

This paper is concerned with the experimental identification of the whole dynamic elastic stiffness tensor of a transversely isotropic clayrock from a single cylindrical sample under loading. Measurement of elastic wave velocities (pulse at 1 MHz), obtained under macroscopically undrained triaxial loading conditions are provided. Further macroscopic (laboratory scale) interpretation of the velocity measurements is performed in terms of (i) dynamic elastic parameters ; and (ii) elastic anisotropy. Experiments were performed on a Callovo-Oxfordian shale, Jurassic in age, recovered from a depth of 613 m in the eastern part of Paris basin in France.

Moreover, a physically-based micromechanical model is developed in order to quantify the damaged state of the shale under loading through macroscopic measurements. This model allows for the identification of the pertinent parameters for a general transversely isotropic orientational distribution of microcracks, superimposed on the intrinsic transverse isotropy of the rock. It is directly inspired from experimental observations and measurements. At this stage, second-and fourthrank tensors α ij and β ijkl are identified as proper dammage parameters. However, they still need to be explicited in terms of micromechanical parameters for the complex case of anisotropy. An illustration of the protocole of this microstructural data recovery is provided in the simpler case of isotropy. This microstructural insight includes cavities geometry, orientation and fluid-content.

Introduction

Callovo-Oxfordian Shale

Clayrocks, and shales in particular, represent approximately two-third of all sedimentary rocks in shallow earth crustal rocks. In oil and gas drilling operations, shales constitute 80 % of all the drilled sections, mainly because they overlie most hydrocarbon bearing reservoirs (immature shales). Furthermore, several countries are considering clayrocks as possible host lithologies for radioactive waste repository, and therefore carrying out research programs to estimate feasibility of such solution. In this trend, the french agency for radioactive waste management, ANDRA, is evaluating the reliability of the Callovo-Oxfordian layer, Jurassic in age, located in the eastern part of France (around Bure), at a depth ranging from 400 m to 700 m. The static mechanical and transport properties of this clayrock have been largely studied for the past few years (complete bibliography in [START_REF] Escoffier | Caractérisation Expérimentale du Comportement Hydromécanique des Argilites de Meuse/Haute-Marne[END_REF]). However, very few elastic wave velocity measurements, all at atmospheric pressure, were performed, and no velocity data at all under mechanical loading are available in the literature for this rock. Indeed, elastic wave velocity measurements are well known as being very sensitive to damage in rocks [START_REF] Schubnel | Dispersion and anisotropy of elastic waves in cracked rocks[END_REF], in particular to evaluate damage in the EDZ (Excavation Damaged Zone). The main outcome of interest here from these research programs is the transversely isotropic nature of Bure clayrock demonstrated experimentally in several labs for several locations and depths in the layer [START_REF] David | Anisotropy of elastic, magnetic and microstructural properties of callovo-oxfordian shale[END_REF]. The main plane of symmetry is quasi-horizontal and corresponds to a very fine layering (bedding).

Shale Lithologies

In general, the monitoring of elastic velocities in shales under mechanical loading are relatively uncommon in view of the abundance of this lithology among shallow earth crustal rocks [START_REF] Stanley | Attenuation anisotropy in shale at elevated confining pressures[END_REF]. This relative rareness is partly due to: (i) the inherent difficulty of transmitting transducer signals through the waterproof pressure chamber in a loading apparatus ; and, on the other hand, (ii) the specificities in preparing and handling shale samples for mechanical testing (chemical sensitivity to water, extremely low permeability...). Most of the dynamic experimental studies reported in the literature on shale samples were performed under hydrostatic loading conditions [START_REF] Johnston | Ultrasonic P and S wave attenuation in dry and saturated rocks under pressure[END_REF][START_REF] Jones | Ultrasonic velocities in cretaceous shales from the williston basin[END_REF][START_REF] Lo | Experimental determination of elastic anisotropy of berea sandstone, chicopee shale, and chelmsford granite[END_REF][START_REF] Johnston | Seismic anisotropy of shales[END_REF][START_REF] Hornby | Experimental laboratory determination of the dynamic elastic properties of wet, drained shales[END_REF][START_REF] Stanley | Attenuation anisotropy in shale at elevated confining pressures[END_REF]. However, elastic wave velocity measurements on shales are reported by [START_REF] Yin | Acoustic Velocity and Attenuation of Rocks: Isotropy, Intrinsic Anisotropy, and Stress Induced Anisotropy[END_REF] un-der triaxial (and polyaxial) loading, and by [START_REF] Podio | Dynamic properties of dry and water-saturated green river shale under stress[END_REF] under uniaxial loading.

Main Goals

The specific experimental setup available in our lab allows for the simultaneous measurement of five different velocities and two directions of strains on the same sample, under triaxial and pore pressure-controlled conditions of loading. This procedure reduces the number of experiments on differently oriented samples usually needed to identify the dynamic properties of a transversely isotropic rock. It also minimizes the errors due to the particular difference between two samples of the supposedly same lithology or same physical state (stress, hydration and saturation histories since recovery).

The main outcomes of this experiment are : (i) identification of the apparent dynamic stiffness tensor of the Callovo-Oxfordian shale from elastic wave velocity measurements ; (ii) assessment of velocity anisotropy, and its evolution under triaxial loading [START_REF] Thomsen | Weak elastic anisotropy[END_REF]. This last step allows for the quantification of the intrinsic and stress-induced anisotropies, leading eventually to an estimation of the microcracks density and distribution evolutions in the shale sample under loading [START_REF] Sayers | Microcracks-induced elastic wave anisotropy of brittle rocks[END_REF][START_REF] Sayers | Stress-dependent seismic anisotropy of shales[END_REF]).

Experimental Procedure

Experimental Setup

The Rock Physics group of the Ecole Normale Supérieure recently acquired a triaxial cell designed for the exploration of the hydromechanical behavior of shallow earth crustal rocks. This apparatus allows for pore pressure, hydrostatic and deviatoric stress to be applied independently on a cylindrical porous rock sample (up to ∅ 40 mm × L 80 mm). The pore, confining and axial pressures may reach 100, 300 and 800 MPa (on the largest samples), respectively. The originality of this loading cell is to allow for a maximum of 32 waterproof signal wires through the wall of the pressure chamber. The control and data acquisition are performed by means of a dedicated software designed in Labview TM environment. The lab is temperature-controlled with an accuracy of ± 0.5 • C.

Sample Description & Preparation

The tested shale was provided by the french agency for radioactive waste management, ANDRA. The received large cores (∅ 90 mm × L 280 mm) were recovered in 1995 from a depth of 613 m in the eastern part of Paris basin, near Bure, in France (MSE 101 borehole). These cores were since preserved in a so-called T1 cell, isolated from gas exchange with the storing environment, and maintained under small, almost hydrostatic, confining stress (. 1 MPa). The sample is dry-cored perpendicularly to the horizontal bedding, ground for parallel end faces, and preserved in waterproof membranes. It is 60 mm long with a 30 mm diameter. For the few days before testing, during the phase of transducers gluing and curing, the sample is maintained in a 100% relative humidity (RH) environment in order to avoid drying and to ensure reproducible initial conditions for different samples with different hydric histories (as much as it is possible with clayrocks). Usually, after putting it in 100% RH atmosphere, the mass of the sample increases then levels off. However, it is assumed that the sample may still be superficially unsaturated. It is tested under macroscopically undrained conditions within the loading cell.

Transducers Arrangement

For the experiment reported here, 22 signal wires were used. Both axial and circumferential strain gages and nine piezoelectric transducers were glued on the cylindrical sample. However, the strain gages signals are rapidly lost at the beginning of this test. The specific arrangement of piezoelectric transducers allows for the measurement of five different elastic velocities as shown in Figure 1. Indeed, let us define the reference frame (x 1 , x 2 , x 3 ), with (x 1 , x 2 ) being the horizontal bedding plane of symmetry of the transversely isotropic cylindrical shale sample. Then, measurements of three compressional velocities and two shear velocities, in three different directions are provided. These velocities are referenced with respect to the bedding plane, i.e., V P (0 • ) for the bedding-parallel compressional velocity, V P (45 • ) for the 45 • -to-bedding compressional velocity, V P (90 • ) for the bedding-perpendicular compressional velocity, V SH (0 • ) for the horizontally polarized, bedding-parallel shear velocity, and V SV (0 • ) for the vertically polarized, bedding-parallel shear velocity. In the case of transverse anisotropy, this latter is equal to the shear wave velocity V S (90 • ) perpendicular to the bedding. Note that the specifically designed compressional piezoelectric transducers are virtually load insensitive. On the other hand, the shear piezoceramics are simply glued to a fitting metal piece, itself glued on the sample, and are therefore load sensitive. Indeed the amplitude of the transmitted signal between two transducers generally increases with external pressure load due to a better contact between a transducer and the sample surface.

Measurement Technique

For velocity measurement, the classical Ultrasonic Pulse Transmission technique is used [START_REF] Birch | The velocity of compressional waves in rocks to 10 kilobars[END_REF][START_REF] Tosaya | Acoustical Properties of Clay-Bearing Rocks[END_REF][START_REF] Yin | Acoustic Velocity and Attenuation of Rocks: Isotropy, Intrinsic Anisotropy, and Stress Induced Anisotropy[END_REF]. This method assumes straight raypaths for the propagation of the elastic perturbation in the medium. It consists in measuring the travel time of a solitary elastic pulse through the rock sample of known travelling wave path length. The frequency of the solitary pulse in this experiment is 1 MHz.

Test Conditions & Loading Path

The chemical sensitivity of this shale to external water and its very low permeability precluded any pore pressure control of the test. Therefore, this test is performed in macroscopically undrained conditions within the triaxial cell (in the hydric state obtained after putting the sample in a 100% RH atmosphere during the sample setting up). The in-situ hydrostatic stress is estimated from depth (' 613 m) and rock density (' 2.45 kg/m 3 ). Indeed, we assume hydrostatic state of stress in this large sedimentary basin, neglecting any deviatoric stress possibly due to tectonic stress build-up. This in-situ stress is estimated around 15 MPa.

The loading starts with a confining pressure cycle between 0 and 55 MPa, i.e., (i) loading from 0 to 55 MPa ; then (ii) unloading from 55 to 15 MPa (overconsolidated test, see Figure 13). The loading rate is 0.013 MPa/s during the pressure-controlled hydrostatic loading. Each point of the experimental curves represents a loading step at which all sample parameters are stabilized before velocity measurements are performed (' 1 hour).

At 15 MPa confining pressure, a deviatoric stress is applied in two steps, i.e., (i) loading from 0 to 46 MPa ; (ii) unloading from 46 to 0 MPa ; then (iii) reloading (Fig. 13). However, the deviatoric loading is displacement-controlled in order to avoid sudden sample failure. Again, velocity measurements at any given loading stage are made after stabilization of all sample parameters. In the specific experiment described here, there are no strain data available during the hydrostatic loading because of experimental difficulties in strain gages gluing. However, an external LVDT allows for the estimation of the sample global axial strain during the deviatoric loading.

Direct Experimental Results

Error Analysis

For this particular experiment, corrections due to the rock deformation in the estimation of the velocity change during the hydrostatic loading are not feasible with accuracy because of a lack of strain measurement. However, a previous experiment on a similar sample, under the same conditions, allowed for strain data acquisition. Therefore, in the hydrostatic part of the present experiment, these strain data are used to correct velocity data. Both samples were cored exactly at the same depth, one centimeter from each other, and no apparent heterogeneity was noted on either of them.

Furthermore, a complete error analysis of the Ultrasonic Pulse Transmission technique has been carried out by [START_REF] Yin | Acoustic Velocity and Attenuation of Rocks: Isotropy, Intrinsic Anisotropy, and Stress Induced Anisotropy[END_REF] or [START_REF] Hornby | Experimental laboratory determination of the dynamic elastic properties of wet, drained shales[END_REF]. This analysis includes the errors in the travel time picks for the measured and reference travel times, and the error in the sample size measurements. In our experimental case, the relative error reaches a maximum of 1.1 % in the velocity measurements, depending on the wave type and propagation direction.

Experimental Data: Elastic Wave Velocities

Figure 2 shows the evolution of the five measured velocities V P (0 • ), V P (45 • ), V P (90 • ), V SH (0 • ) and V SV (0 • ), as a function of confining pressure varying between 0 and 55 MPa. All velocities increase with increasing pressure, and a hysteresis effect exists upon pressure release from 55 MPa toward 15 MPa. Note that in the 0 -55 MPa pressure range, the velocity variations due to pressure (subscript p) are ordered, i.e., ∆ p V P (0

• ) < ∆ p V P (45 • ) < ∆ p V P (90 • ).
Even the amplitudes of the hysteresis velocity lags (subscript h) are ordered in the same way, i.e., ∆ h V P (0

• ) < ∆ h V P (45 • ) < ∆ h V P (90 • ). Measured initial ve- locities at atmospheric pressure are V o P (0 • ) = 3120 m/s, V o P (45 • ) = 2771 m/s, V o P (90 • ) = 2314 m/s, V o SH (0 • ) = 1906 m/s and V o SV (0 • ) = 1693 m/s.
Figure 3 shows the axial stress-axial strain curve during the deviatoric part of the loading path. Two loading cycles have been performed, and after the first unloading, irreversible strain is observed. The axial loading (displacementcontrolled) reaches a peak stress σ peak = 50 MPa at 1.9% axial strain, before decreasing toward σ end = 42 MPa and stabilizing at this value. Post-mortem analysis of the sample showed a localization of the deformation in a shear band inclined at an angle of θ ' 45 • with respect to the horizontal bedding plane.

Figure 4 shows the evolution of the five measured velocities

V P (0 • ), V P (45 • ), V P (90 • ), V SH (0 • ) and V SV (0 •
), as a function of increasing axial load (at 15 MPa confining pressure), applied in two steps, until peak stress σ peak = 50 MPa. Note the monotonic decrease of V P (0 • ) and V SH (0 • ). Conversely, V P (90 • ) is monotonically increasing until peak stress. However, V P 45 (45 • ) and V SV (0 • ) first increase then noticeably decrease after 30 MPa axial stress. Moreover, for all measured velocities, a residual change in velocity is observed upon axial unloading from 46 toward 0 MPa, which typically denotes an irreversible change in the (micro)structure of the rock. Further reloading of the sample produces qualitatively the same velocity evolutions as in the first loading step, until maximum axial stress is reached (50 MPa).

4 Macroscopic Theoretical Interpretation

Shale Qualitative Description

Clayrocks, and shales in particular, are usually considered as a multi-scale medium composed of very anisotropic clay platelets, surrounding possible inclusions of other stiffer minerals (mica, pyrite) or kerogen. These clay platelets are grouped together to form so-called clay particles, which form a more or less complex network, defining micro-, meso-and macro-scale Representative Elementary Volumes (R.E.V.) (Fig. 5). The more the clay particles are aligned at the mesoscale, the more the clayrock is anisotropic at the macroscale. In particular, our Jurassic shale is recognized to be transversely isotropic at the macroscopic scale [START_REF] David | Anisotropy of elastic, magnetic and microstructural properties of callovo-oxfordian shale[END_REF].

Dynamic Stiffness Matrix

Definition

The previous qualitative description allows for the modeling of this rock as transversely isotropic (hexagonal symmetry) at the macroscopic scale, with x 3 as the symmetry axis, and the plane (x 1 , x 2 ) as the plane of isotropy. This type of materials may be adequately described in the elastic regime by five independent stiffness constants relating stress to strain components. They form a fourth order tensor which may be written in the Voigt notation as a twelve components matrix in the reference frame (x 1 , x 2 , x 3 ), i.e.,

                 σ 11 σ 22 σ 33 σ 23 σ 13 σ 12                  =                  C 11 C 12 C 13 0 0 0 C 12 C 11 C 13 0 0 0 C 13 C 13 C 33 0 0 0 0 0 0 C 44 0 0 0 0 0 0 C 44 0 0 0 0 0 0 C 66                                   ε 11 ε 22 ε 33 2ε 23 2ε 13 2ε 12                  , (1)
where C 12 = C 11 -2C 66 . Note that here σ ij and ε ij stand for the stress and strain average tensors at the laboratory (macroscopic) scale (Fig. 5).

When evaluated from elastic wave velocity measurements, this tensor characterizes the dynamic apparent behavior of the material. These five independent C ij constants may be related uniquely to the five classical engineering elastic constants of a transversely isotropic medium, E 1 , E 3 , ν 12 , ν 13 and ν 31 for instance [START_REF] Podio | Dynamic properties of dry and water-saturated green river shale under stress[END_REF]), i.e.,

E 1 = D C 11 C 33 -C 2 13 , E 3 = D C 2 11 -C 2 12 , ν 12 = C 12 C 33 -C 2 13 C 11 C 33 -C 2 13 , ν 13 = C 13 (C 11 -C 12 ) C 11 C 33 -C 2 13 , (2)
ν 31 = C 13 C 11 + C 12 .
where D is the following determinant

D = ¯C11 C 12 C 13 C 12 C 11 C 13 C 13 C 13 C 33 ¯.
(3)

These five independent C ij constants are related to the five elastic wave velocities measured experimentally as follows (see Figure 1 for geometrical reference configuration):

C 11 = ρV 2 P (0 • ) , C 33 = ρV 2 P (90 • ) , C 44 = ρV 2 SV (0 • ) , C 66 = ρV 2 SH (0 • ) , (4) C 13 = -C 44 + [(C 11 + C 44 -2ρV 2 P (45 • )) (C 33 + C 44 -2ρV 2 P (45 • ))]

Dimensionless Anisotropy Parameters

Definition

In the specific context of anisotropy analysis at the laboratory scale, it is convenient to express three of the five dynamic elastic constants of (1) in terms of three anisotropy parameters, initially introduced by Thomsen (1986), i.e.,

ε = C 11 -C 33 2 C 33 , γ= C 66 -C 44 2 C 44 , δ = (C 13 + C 44 ) 2 -(C 33 -C 44 ) 2 2 C 33 (C 33 -C 44 ) , (5)
where ε measures the P-wave anisotropy, γ measures the S-wave anisotropy and δ is totally independent of horizontal velocity and may be either positive or negative for geomaterials. Actually, δ quantifies the dynamic effect of the fluid phase contained in the macroscopic R.E.V., i.e., pore water and/or air. An additional anellipticity parameter may be defined as ε-δ, which quantifies the departure of the slowness surface from an ellipsoidal shape [START_REF] Thomsen | Weak elastic anisotropy[END_REF][START_REF] Hornby | Experimental laboratory determination of the dynamic elastic properties of wet, drained shales[END_REF][START_REF] Schubnel | Dispersion and anisotropy of elastic waves in cracked rocks[END_REF].

The main advantage of these parameters is their dimensionless nature, therefore facilitating qualitative insight into the material behavior.

Evolution of Anisotropy under Loading

Figures 10 and 11 show the evolution of velocity anisotropy parameters ε, γ and δ as a function of the hydrostatic and deviatoric loadings, respectively.

Note the significant decrease of all anisotropy parameters with increasing hydrostatic loading. All anisotropy parameters seem to level-off, with non-zero values at the highest confining pressure, i.e.,

ε i = 0.21, γ i = 0.07, δ i = 0.28. (6)
Note that this so-called "intrinsic" (residual) anisotropy at high confining pressure represents half of the total anisotropy at ambient pressure (Fig. 10).

Moreover, ε and δ seem to be very sensitive to confining pressure, however, under axial stress, δ seems to be the most sensitive parameter. Anisotropy globally reduces under axial load.

5 Diffuse Damage Modeling

Main Assumptions

At the microscale (Fig. 5), in the particular case of clays, the possible damage may be interpreted in terms of contact regions between clay platelets (Fig. 12).

In the following developments, the words crack or microcrack are to be understood as crack-like voids which may be close to an ellipsoidal shape of characteristic half lengths, a and w, with w ¿ a. Indeed, "crack-like" has to be understood as "strongly oblate". In that case, it may be shown using Hill's auxiliary theorem that further details on the exact geometry of the void are irrelevant, as far as the elastic properties are concerned [START_REF] Hill | Elastic properties of reinforced solids: Some theoretical principles[END_REF][START_REF] Sevostianov | Recovery of information on the microstructure of Porous/Microcracked materials from the effective Elastic/Conductive properties[END_REF].

Damage is strongly dependent on the stress state of the rock. As shown by the anisotropy evolution curves (Fig. 10), at high enough confining pressure, anisotropic voids may close while anisotropy reduces. However, a residual anisotropy remains, characterizing an intrinsic textural/mineral anisotropy.

Moreover, after the initial recovery of the shale core from the borehole, microcracks may appear due to stress release (from in-situ stress to atmospheric pressure). These microcracks tend to be aligned within the bedding plane because of the larger compliance of the rock in the direction perpendicular to this. According to classical rock mechanics theory, further deviatoric loading induces microcracking parallel to the major principal stress (parallel to the axial load). Figure 13 summarizes these remarks and illustrates the transverse isotropy of microcracks orientational distribution under either hydrostatic or deviatoric loading.

Eventually, due to the very low permeability of this shale and the high domi-nant frequency of the induced elastic wave (f exp = 1 MHz), fluid pressure relaxation around the travelling wave path may not have time to occur. In order to quantify the extent of this assumption, we borrow arguments from dynamic poroelasticity theory. Indeed [START_REF] Biot | Generalized theory of acoustic propagation in porous dissipative media[END_REF] introduced, on a phenomenological basis, the limit Biot frequency f B at which inertial and viscous fluid effects are of the same magnitude, i.e.,

f B = ηφ 2πkρ f , (7)
where

ρ f h FT 2 /L 4 i is the fluid density, k [L 2
] denotes the macroscopic average permeability of the porous medium, η [FT/L 2 ] is fluid dynamic viscosity, φ stands for the porosity.

If f exp ¿ f B , then fluid is mostly in the viscous regime and inertial effects are negligible at the mesoscale (Fig. 5). In particular, by means of published experimental data for our shale, the Biot frequency is estimated to be f B ' 2 10 2 to 2 10 4 [GHz] in the case of water saturating fluid (reported data for this shale are: k ' 1 to 100 [10 -21 m 2 ], φ ' 10%). This frequency f B is much larger than the experimental perturbation frequency, which allows us to assume that water saturating fluid is in the viscous regime when the microscopic R.E.V. is crossed by the travelling elastic wave. This approach is clearly a rough approximation since we did not provide any information about the geometry of the void space in the poroelastic model (pores are treated as isotropic by default at the mesoscale R.E.V.). Indeed, the presence of cracklike pores (micro-anisotropy), or gas partial saturation (micro-heterogeneity) may modify strongly this first conclusion. Furthermore, these causes may also induce squirt-flow extra dissipation phenomena, both in the saturated and unsaturated cases. However, this approach illustrates the complex nature of the behavior of such porous materials. Therefore, we need to build a robust model describing this medium in a specifically particular condition, from which we may evaluate the deviation of experimental data from this physically wellbased specific state. Indeed, we will focus, in a first modeling step, on the case of high frequency elastic perturbation. In other words, whatever the fluid contained in anisotropic crack-like voids, it is microscopically a "frozen fluid" in ellipsoidal voids. This type of model allows for the estimation of the density and orientational distribution of microcracks in the medium (diffuse damage) at a given stress state. Since elastic wave velocities are known to be sensitive to discontinuities, it seems natural to evaluate these model parameters by means of this experimental technique [START_REF] Kachanov | Elastic solids with many cracks and related problems[END_REF][START_REF] Sayers | Microcracks-induced elastic wave anisotropy of brittle rocks[END_REF][START_REF] Sayers | Stress-dependent seismic anisotropy of shales[END_REF].

Analytical Effective Medium Modeling

Under the above assumptions, the model takes into account the intrinsic transverse isotropy of the rock matrix (without strongly anisotropic voids). A transversely isotropic orientational distribution of these strongly anisotropic voids is then superimposed. This last step implicitly requires the non interacting condition between different voids in order to end up with an analytical model [START_REF] Kachanov | Elastic solids with many cracks and related problems[END_REF][START_REF] Sayers | Microcracks-induced elastic wave anisotropy of brittle rocks[END_REF][START_REF] Sayers | Stress-dependent seismic anisotropy of shales[END_REF]. This voids-induced transverse isotropy has the same principal directions as the intrinsic symmetry under either hydrostatic or deviatoric loading.

We use a compliance formulation of elasticity equations in order to discriminate between the two sources of anisotropy, i.e., intrinsic and voids-induced anisotropies

ε ij = ³ S o ijkl + ∆S ijkl ´σkl , (8)
where S o denotes the intrinsic compliance tensor of the shale (clay minerals transverse isotropy) ; ∆S is the damage-induced excess compliance tensor ; ε and σ are the average strain and stress second-rank tensors at the microscale, respectively.

Assuming a condition of non-interacting voids, for each anisotropic void the average displacement discontinuity vector [u] may be related to the average stress field σ through the compliance tensor B as follows

[u i ] = B ij σ jk n k , ( 9
)
where n is the unit vector perpendicular to the discontinuity (Fig. 12).

As pointed out by [START_REF] Sayers | Microcracks-induced elastic wave anisotropy of brittle rocks[END_REF], the compliance tensor of one discontinuity in a general infinite elastic medium may be expressed in terms of its spherical and deviatoric components B N and B T , respectively, as follows

B ij = B N n i n j + B T (δ ij -n i n j ) , (10)
where δ ij stands for the Kronecker delta tensor. Therefore, following the rationale developed by [START_REF] Sayers | Microcracks-induced elastic wave anisotropy of brittle rocks[END_REF], the extra compliance tensor ∆S due to many non-interacting crack-like voids at the macroscale (all voids subjected to the same stress field σ) is

∆S ijkl = 1 V X r B r T 4 A r ³ δ ik n r j n r l + δ il n r j n r k + δ jk n r i n r l + δ jl n r i n r k + 1 V X r (B r N -B r T ) A r n r i n r j n r k n r l , (11)
where superscript r stands for the r th contact, and A r denotes the r th crack surface area. Defining a second-and fourth-rank tensors α and β as

α ij = 1 V X r B r T A r n r i n r j , (12a)
β ijkl = 1 V X r (B r N -B r T ) A r n r i n r j n r k n r l , (12b)
leads to the following expression of the extra compliance tensor at the macroscale

∆S ijkl = 1 4 (δ ik α jl + δ il α jk + δ jk α il + δ jl α ik ) + β ijkl . ( 13
)
Note that tensors α and β are symmetric with respect to all rearrangements of the indices. For a general transversely isotropic distribution of discontinuities, α 11 = α 22 , β 1111 = β 2222 , and β 1212 = β 1122 = β 1111 /3. Therefore, (13) yields

∆S 11 = ∆S 22 = α 11 + β 1111 , ( 14
) ∆S 33 = α 33 + β 3333 , (15)
∆S 66 = 2α 11 + 4 3 β 1111 , (16)
∆S 44 = ∆S 55 = α 11 + α 33 + β 1133 , ( 17
) ∆S 13 = ∆S 23 = β 1133 , ( 18
)
in Voigt double index convention for the tensor ∆S.

Damage Parameters Identification

First, we estimate the intrinsic dynamic stiffness tensor of the shale from velocity measurements at the highest confining pressure available (55 MPa), i.e., when the crack-like voids are assumed to be closed. This state is taken to be a reference at which all crack-like voids are nominally closed. The existence of these voids in the sample is induced by the hydromechanical history following recovery from deep drilling operation: (i) in-situ porosity, (ii) porosity induced by stress relief ; and (iii) pore fluid overpressures during and after recovery. Therefore, any variation of the apparent dynamic stiffnesses at confining pressures lower than 55 MPa (reference state) is supposed to be due to voids geometry evolution (opening/closing of porosity). This extra stiffness tensor is inverted to yield the extra compliances and the system ( 14)-( 18) is solved for α 11 , α 33 , β 1111 , β 3333 and β 1133 at all loading stages, which yields

α 11 = -2∆S 11 + 3 2 ∆S 66 , (19)
α 33 = 2∆S 11 -4∆S 13 + ∆S 44 - 3 2 ∆S 66 , (20)
β 1111 = 3∆S 11 - 3 2 ∆S 66 , (21)
β 3333 = -2∆S 11 + 4∆S 13 + ∆S 33 -∆S 44 3 2 ∆S 66 , (22)
β 1133 = ∆S 13 . ( 23
)
Figures 14 and15 show the evolution of the damage parameters as a function of the hydrostatic and deviatoric loadings, respectively. Under increasing confining pressure, as expected, there is a global decrease in the absolute value of most damage parameters toward the nominal zero value at 55 MPa, even if α 33 is always larger than α 11 . Indeed, at atmospheric pressure (after initial recovery), there are more horizontally oriented crack-like voids than vertically oriented, and under hydrostatic load, both families do close (Fig. 14). However, note the increase in the density of crack-like voids perpendicular to bedding α 11 under deviatoric loading at 15 MPa confining pressure, while α 33 seems stable until 45 MPa axial stress, then decreases drastically between 45 and 50 MPa, illustrating the further closure of bedding-parallel cracks (Fig. 15).

The precise analysis and interpretation of the evolutions of β 1111 , β 3333 and β 1133 is still under investigation in the anisotropic case. They are related to the geometry, orientation and fluid-content of the crack-like voids.

These micromechanics-based crack parameters are directly related to macroscopically measurable (phenomenological) quantities. More precisely, α ij are the parameters sensitive to discontinuities, and β ijkl take into account the properties of the fluid possibly contained within the pore volume (air, water or both). However, they do not yet reflect an explicit link with the physical properties of this filling fluid at the mesoscale (Fig. 5). Indeed, they are expressed in terms of B N and B T , which are lump parameters including the solid and the fluid phase elastic properties, as well as the geometry of the discontinuity (equant and/or crack-like pores) within the mesoscale R.E.V.. The explicit dependence of B N and B T on these mesoscale quantities (fluid and solid properties) is still under investigation in the case of an intrinsically transversely isotropic medium with arbitrary distribution of crack-like voids.

However, for illustration purposes, let us explicit the link between B N and B T and the microscopic physical parameters in the simpler and well-documented case of an isotropic elastic solid with randomly oriented and located fluid-filled crack-like pores. Indeed, in this case, the components B N and B T of (10) for a single crack are [START_REF] Shafiro | Materials with fluid-filled pores of various shapes: Effective elastic properties and fluid pressure polarization[END_REF][START_REF] Sevostianov | Recovery of information on the microstructure of Porous/Microcracked materials from the effective Elastic/Conductive properties[END_REF])

B N = 16a (1 -ν 2 ) 3πE µ 1 - 1 1 + λ ¶ , (24a)
B T = 32a (1 -ν 2 ) 3π (2 -ν) E , (24b)
where E and ν denote Young's modulus and Poisson's ratio of the elastic solid, respectively ; a stands for the radius of the penny-shaped void ; and λ is a solid-fluid coupling term explicited below. It is a function of the filling fluid elastic constant K f [MPa] and the aspect ratio of the void ζ = w/a. The explicit relation is [START_REF] Shafiro | Materials with fluid-filled pores of various shapes: Effective elastic properties and fluid pressure polarization[END_REF][START_REF] Sevostianov | Recovery of information on the microstructure of Porous/Microcracked materials from the effective Elastic/Conductive properties[END_REF])

λ = ζπ 4 (1 -ν 2 ) Ã E K f -3 (1 -2ν) ! , (25)
where w stands for the half-opening of the crack-like void.

Note that in the case of a dry cavity, the second term in the brackets of (24a) cancels (λ → ∞ while the incompressibility of the fluid K f → 0), and B N no longer depends on the aspect ratio of the cavity ζ.

For a set of N non-interacting, identical, randomly oriented and located, fluid-filled cracks in an isotropic elastic solid, α ij and β ijkl are isotropic and symmetric with respect to all rearrangements of indices. Therefore, (12a) and (12b) become

α ij = πNa 2 3V B T δ ij , (26a)
β ijkl = πNa 2 15V (B N -B T ) (δ ij δ kl + 2 (δ ik δ jl + δ il δ jk )) . (26b)
where we recall that δ ij stands for the Kronecker delta tensor.

And the non-vanishing components of α ij and β ijkl in this simple isotropic case are

α 11 = α 22 = α 33 = 1 3 α kk , β 1111 = β 2222 = β 3333 , (27)
β 1122 = β 1133 = β 2233 = β 1313 = β 2323 = 1 3 β 1111 .
Therefore,

α 11 = πN a 2 3V B T , (28a)
β 1111 = πN a 2 3V (B N -B T ) . (28b)
By inserting (24a) and ( 24b) in ( 28a) and (28b), we notice that the model is mainly function of two independent damage parameters, i.e., ρ and ρ 1 . Among these two parameters, we recognize the classical scalar crack density parameter for identical cracks

ρ ρ = 1 V N X r (a r ) 3 = Na 3 , (29)
and on the other hand, we have an additional parameter ρ 1 depending on the fluid content (fluid elstic properties and saturation of the porous volume)

ρ 1 = 1 V N X r (a r ) 3 1 + λ r = N a 3 1 + λ . (30)
In the more general case of transverse isotropy, further studies are needed in order to reach such a clear insight on the physical problem. At this stage, only "lump" damage parameters α ij and β ijkl have been inverted from experimental data in this more complex anisotropic problem, however, the illustrated example in the simpler case of isotropy shows the protocole in recovering microstructural data from the macroscopic measurements.

Discussion

Anisotropy and Damage Issues

The main factors contributing to velocity anisotropy of rocks are: (i) preferred orientation of void space (pores and cracks) ; (ii) textural-structural features (bedding or foliation) ; and (iii) anisotropic constituent minerals. Damageinduced anisotropy reduces with increasing hydrostatic load, leading eventually to the estimation of the intrinsic anisotropy at high confining pressure, due to factors (ii) and (iii) only. If we assume that crack-induced and intrinsic anisotropies are additive [START_REF] Vernik | Velocity anisotropy in shales: A petrophysical study[END_REF], they may be discriminated from each other in the experimental data.

The global analysis of the anisotropy data leads to a conceptual model of the damage history of a shale after coring and recovery. In general, (i) stress relief, (ii) micro-hydraulic fracturing, (iii) bottomhole stressing, (iv) dessication, and (v) in-situ microcracks, contribute to the final state of microcracking in a rock sample in the lab. It is likely that after initial core recovery, microcracks will mostly open in a plane aligned with the bedding plane, perpendicularly to the direction of larger compliance of the transversely isotropic rock. This fact is experimentally verified for our Jurassic shale. Indeed, microcracks close under confining pressure, preferentially in the direction perpendicular to the bedding (bedding-parallel microcracks), since α 33 > α 11 . The deviatoric loading induces (i) further closure of bedding-parallel microcracks, starting from a certain level of axial stress from which α 33 still decreases drastically ; and (ii) additional damage (bedding-perpendicular microcracks), since α 11 increases with increasing load. The non-monotonic variation of α 33 between 30 and 40 MPa axial stress (Fig. 15) need to be confirmed by further experiments. However, if this is confirmed, this effect may be due to the presence of liquid and variation of saturation state under load.

Fluid Content, Saturation and Effective Stress Issues

As pointed out by [START_REF] Jones | Ultrasonic velocities in cretaceous shales from the williston basin[END_REF], the partial saturation of the sample manifests itself by the velocity hysteresis under cyclic confining pressure. They noticed that this hysteresis is drastically reduced when the shale sample is re-humidified, and even cancels at very high confining pressure (400 MPa). Indeed, at ambient pressure, the void space is partially gas-saturated, in particular at the boundaries of the sample where gas exchange with atmosphere may occur (surface drying). When submitted to a hydrostatic load, the microcracks closure reduces the void space, which is eventually totally occupied by the pore water at very high confining pressure. In our experiment, it seems that even at 55 MPa confining pressure, elastic wave velocities hysteresis upon pressure release still happen (Fig. 2). This fact suggests that even at 55 MPa confining pressure, the pore space is still not fully saturated by water.

This brings out the question of effective pressure in the sample under loading. We will assume that during the hydrostatic loading, the confining pressure is equal to the poroelastic effective pressure, since the pore space is not water-saturated and the incompressibility of gas is negligible compared to that of the liquid. The sample acts as under drained conditions. However, identification of the solid-fluid coupling parameter δ from experimental data may lead to a quantification of these changes in saturation along the test.

Concluding Remarks

Anisotropic elastic properties of a quasi-saturated Jurassic shale under macroscopically undrained conditions were determined by means of ultrasonic ve-locity measurements. These five velocity measurements were performed on the sample while under hydrostatic and deviatoric loadings in a triaxial cell. The experimental method is based on a single core test, minimizing the errors due to the particular difference between two samples of the supposedly same lithology or same physical state (stress and saturation history since recovery). However, in the present study, strain data during the hydrostatic phase of the mechanical loading have been obtained on a different sample cored very close to the sample of interest. Indeed, the strain gages glued on the sample of interest did not work properly at the time of the main experiment.

Analysis and quantification of anisotropy allowed for the estimation of pertinent parameters related to damage evolution. Further development of the model may lead to an estimation of either fluid content or void average aspect ratio, or even both.

The first practical application of this type of experiment and micromechanical modeling is to provide input parameters for modeling wave propagation in porous, fluid-infiltrated, very low permeability shales. A natural phenomenological way to model such medium is to use Biot's dynamic theory of poroelasticity, possibly extended to account for squirt-flow phenomenon due to the existence of an equant porosity connected to an anisotropic one (crack-like shape).

Acknowledgements

Support and shale samples for this research program have been provided by the french agency for radioactive waste management, ANDRA. This support is gratefully acknowledged. I would also like to thank Guy Marolleau and Thierry Descamps for their crucial help in the setting of the experiment in the triaxial cell.

University, Stanford, USA.

List of Figures

Mean Stress

Differential Stress

Y i e l d C u r v e D i l a t a n c y

In-situ State

Over-Consolidated State

Ambiant Conditions (after initial recovery)

Fig. 2 .

 2 Fig. 2. Evolution of compressive and shear wave velocities in the Callovo-Oxfordian shale under hydrostatic loading between 0 and 55 MPa.

Fig. 3 .Fig. 4 .Fig. 5 .

 345 Fig. 3. Axial stress-axial strain curve at 15 MPa confining.

Fig. 6 .

 6 Fig. 6. Evolution of dynamic stiffness parameters of the Callovo-Oxfordian shale under hydrostatic loading between 0 and 55 MPa.

Fig. 7 .Fig. 8 .

 78 Fig. 7. Evolution of dynamic stiffness parameters of the Callovo-Oxfordian shale under deviatoric loading between 0 and 50 MPa, at 15 MPa confining pressure.

Fig. 9 .Fig. 11 .

 911 Fig. 9. Evolution of dynamic engineering constants of the Callovo-Oxfordian shale under deviatoric loading between 0 and 50 MPa, at 15 MPa confining pressure.

  Fig.12. A region within a shale showing a local alignement of clay platelets (e.g., clay particle). Contacts between clay platelets or between clay particles may be modeled as crack-like discontinuities (after[START_REF] Sayers | Stress-dependent seismic anisotropy of shales[END_REF]).

Fig. 13 .Fig. 14 .Fig. 15 .

 131415 Fig. 13. Conceptual model of the evolution of anisotropic voids in a transversely isotropic shale (intrinsic and voids-induced anisotropies) under general loading.

  1/2 .4.2.2 Evolution of Shale Stiffness under LoadingFigures 6 and 7 illustrate the evolution of the dynamic stiffness parameters C ij as a function of hydrostatic and deviatoric loadings, respectively. Measured initial dynamic moduli at atmospheric pressure are C o 11 = 23.8 GPa, C o 33 = 13.1 GPa, C o 44 = 7 GPa, C o 66 = 8.9 GPa and C o 13 = 3.8 GPa. Note that Figure 6 allows for the estimation of the crack-free dynamic moduli at low confining pressures by a linear extrapolation of the high pressure moduli.Equivalently, figures 8 and 9 illustrate the evolution of the dynamic engineering constants, i.e., Young's moduli E 1 and E 3 , and Poisson's ratio ν 12 , ν 13 and ν 31 , as a function of hydrostatic and deviatoric loadings, respectively.

	Measured initial dynamic engineering constants at atmospheric pressure are
	E o 1 = 21.7 GPa, E o 3 = 12 GPa, ν o 12 = 0.22 , ν o 13 = 0.23 and ν o 31 = 0.13.

1

  Schematics of the positions of the piezoceramic transducers with respect to the sample (with horizontal bedding plane).

	12 A region within a shale showing a local alignement of clay
		platelets (e.g., clay particle). Contacts between clay platelets or between clay particles may be modeled as crack-like P-transducer
	This figure also illustrates the polarization of the shear discontinuities (after Sayers (1999)). 13 Conceptual model of the evolution of anisotropic voids in Polarized S-transducer	30
	transducers, as well as the deduced elastic wave velocities. a transversely isotropic shale (intrinsic and voids-induced and 55 MPa. shale under hydrostatic loading between 0 and 55 MPa. Callovo-Oxfordian shale under hydrostatic loading between 0 14 Evolution of damage parameters in the Callovo-Oxfordian 2 Evolution of compressive and shear wave velocities in the anisotropies) under general loading. Elastic Wave Propagation	23 23 32 31
	3 15 Evolution of damage parameters in the Callovo-Oxfordian Axial stress-axial strain curve at 15 MPa confining. shale under deviatoric loading between 0 and 50 MPa, at 15	24
	4	Evolution of compressive and shear wave velocities in the MPa confining pressure.	3 3
		Callovo-Oxfordian shale under deviatoric loading between 0
		and 50 MPa, at 15 MPa confining pressure.	2 4
	5	Multiscale R.E.V.'s illustrating the geometric arrangement of a
		clay matrix. Note that Callovo-Oxfordian shale is transversely
		isotropic at the macroscopic scale.	25
	6	Evolution of dynamic stiffness parameters of the Callovo-
		Oxfordian shale under hydrostatic loading between 0 and 55
		MPa.	25
	7	Evolution of dynamic stiffness parameters of the Callovo-
		Oxfordian shale under deviatoric loading between 0 and 50
		MPa, at 15 MPa confining pressure.	2 6
	8	Evolution of dynamic engineering constants of the Callovo-
		Oxfordian shale under hydrostatic loading between 0 and 55
		MPa.	27
	9	Evolution of dynamic engineering constants of the Callovo-
		Oxfordian shale under deviatoric loading between 0 and 50
		MPa, at 15 MPa confining pressure.	2 8
	10 Evolution of elastic anisotropy parameters of the Callovo-
		Oxfordian shale under hydrostatic loading between 0 and 55
		MPa.	29
	11 Evolution of elastic anisotropy parameters of the Callovo-
		Oxfordian shale under deviatoric loading between 0 and 50
		MPa, at 15 MPa confining pressure.	3 0