Rachid Ababou

Hassane Fatmi

Jean-Michel Matray

Christophe Nussbaum

David Bailly

Statistical Analyses of Pore Pressure Signals in Claystone During Excavation Works at the Mont Terri Underground Research Laboratory

Introduction

In many countries (such as Belgium, Germany, France, Japan, Switzerland, and United Kingdom), deep argillaceous formations are considered as potential host rocks for geological disposal of high-level and intermediate-level long-lived radioactive wastes. Some of these countries are investigating the suitability of high compacted clay-rich rocks at depths down to around 500 m below the ground surface. The general disposal concept comprises a network of drifts and tunnels linked to the surface by shafts and ramps, all artificially ventilated. Research is ongoing in Underground Research Laboratories, like the Mont Terri site in the Swiss Jura, to assess and ensure the safety of the repositories for the full decay life of the radioactive waste, i.e. the capacity of the hypothetical repository to prevent the migration of radionuclides towards the biosphere.

One of the important safety considerations regarding nuclear waste disposal repositories are damages and disturbances induced by excavation works to the rock-mass in the vicinity of the openings. These excavations can alter the initial and favorable isolation properties of the host rock. The most obvious problem is that of changes in the stress field around the tunnel during excavation that induces a fracture network, namely the Excavation Damaged Zone or EDZ, consisting mainly in unloading joints and shrinkage cracks [START_REF] Bossart | Geological and hydraulic characterisation of the excavation disturbed zone in the Opalinus Clay of the Mont Terri Rock Laboratory[END_REF] and therefore characterized by irreversible changes. This phenomenon creates an increase in permeability of potentially several orders of magnitude higher than that of the undamaged host rock [START_REF] Blümling | The excavation damaged zone in clay formations time-dependent behaviour and influence on performance assessment[END_REF]. The EDZ and its raised permeability could possibly affect the performance and safety of the repository as it could act as a preferential pathway, allowing radionuclides to migrate to the surface, bypassing bentonite seals along the drifts.

However, excavation works may also contribute to reversibly disturb the geological medium in the so-called Excavation disturbed Zone (EdZ). This term (EdZ) designates a zone that is thought to be disturbed by elastic hydro-mechanical variations of pressure and stress, but without irreversible changes in flow and transport properties. Within the EdZ, there are no negative effects on the long-term safety, but short-term effects may exist as the response to hydro-mechanical coupling. For instance, mine-by test experiments performed during excavation works often show an increase in pore pressures in the vicinity of the newly excavated gallery. This increase may exceed 1 MPa (10 bars), for instance during excavation of Ga03 in the Tournemire URL [START_REF] Massmann | Modelling of Excavation Induced Coupled Hydraulic-Mechanical Processes in Claystone[END_REF]. This pore pressure increase appears either instantaneous or delayed depending on the position of the measurement sections with respect to the excavation front.

In short, the objective of this text is to study the isolation properties of a geologic repository against the migration of radionuclides through the host rock and eventually towards the biosphere. In this context, the hydro-geologic properties of the host rock must be characterized at different scales (e.g., near field vs. far field).

The permeability is a very important property in terms of isolation, but in this work, we focus on novel methods for identifying two other related properties: (i) specific storativity (as the elastic response to earth tides), and (ii) the effective dynamic porosity (using the barometric effect). Indeed, the analysis of pore pressure and atmospheric pressure time series makes it possible in some cases to determine hydraulic parameters such as specific storage coefficient and effective porosity using simplified groundwater flow and compressibility models (e.g., [START_REF] Marsaud | Estimation des caractéristiques physiques d'aquifères profonds à partir de l'incidence barométrique et des marées terrestres[END_REF]). These models have already been applied to obtain a preliminary characterization of the compacted clay rock of the Tournemire URL in Aveyron, France (Fatmi et al. 2005).

More precisely, our main goal in this text is to characterize the above mentioned hydraulic and mechanical properties of the clay rock from analyses of measured atmospheric and pore pressure signals, and also, to study the response of the clay rock's "Excavation disturbed Zone" ("EdZ") during the excavation of galleries.

The data analysis techniques rely strongly on statistical concepts (theory of random processes) and on the decomposition of signals on various bases:

•

Harmonic basis (sine, cosine): Fourier spectrum S(f) vs. frequency (f).

•

Wavelet multiresolution bases (Daubechies): wavelet components Cj(t) at different time scales "j".

These techniques are used for interpreting the proposed hydromechanical models, in order to identify the hydrogeologic properties of the claystone from measured signals (possibly disturbed by excavation works):

• Wavelet decomposition of pore pressure signals (selection of semi-diurnal earth tide component), and statistical envelope analysis (effect of moving excavation front): earth tide effect model → specific storativity; • Spectral cross-analysis (spectral gain) of air pressure vs. pore pressure signal: barometric effect model → porosity.

The study presented in this text is entirely based on data and time series collected at the excavated geologic claystone site of the Mont Terri URL in the Swiss Jura (Mont Terri Projectinternational consortium), and we focus on claystone properties at "near field" scales, within meters or tens of meters from excavated galleries at the Mont Terri site. Similar techniques are currently being developed for other datasets collected at the Tournemire URL (Aveyron, France) operated by the "Institut de Radioprotection et de Sûreté Nucléaire" (France).

In this work, in order to achieve the above mentioned objectives, we use pressure signals obtained over a period of ten years at Mont Terri (1996Terri ( /2005)), from which records as long as one year or more will be analysed in this text (see Section 4). But, to analyse the effects of a moving excavation front, we also focus on a period of 5 months of "syn-excavation" period for gallery "Ga98", and on a much shorter period of fast pressure changes lasting only a couple of weeks (see Section 5). In total, we study both long and short time scale records, from years to weeks.

Some results of pressure signal pre-processing and long time scale analyses were presented in Fatmi et al. (2007[START_REF] Fatmi | Statistical pre-processing and analyses of hydrogeo-meteorological time series in a geologic clay site (methodology and first results for Mont Terri's PP experiment)[END_REF] at an early stage of this study. The present text provides a significant update of the results of long time scale analyses, and presents also new types of results on the analysis of excavation effects on pore pressures at very short time scales (weeks).

Accordingly, the outline of this chapter is as follows.

Section 2 describes the hydro-geologic site and the pore pressure measurements at the Mont Terri site (in the framework of experiment LP14).

Section 3 is devoted to the mathematical and statistical methods of signal pre-processing (e.g. reconstruction of missing data) and of signal analysis (time-lag correlation functions, frequency spectra and cross-spectra, wavelets, envelopes). For hydro-geologic interpretation of pressure signals, the methods involve, not only statistical treatments, but also simplified elastic hydro-mechanical models: the concept of a specific storativity coefficient "Ss" (in relation to earth tides); and the barometric efficiency (involving both "Ss" and porosity "Φ") in relation to air pressure fluctuations. These models are briefly described at the end of the section (Box N°2) and their application is developed in the remaining sections of the text (see below).

Section 4 develops the implementation of "long time scale" statistical analyses and hydrogeologic interpretations of signals, in the case where the objective is to obtain a global estimate of the two hydro-geologic coefficients (Ss,Φ) in two steps, first "Ss" from semidiurnal earth tide effects using multi-resolution wavelet analysis, and secondly "Φ" using the barometric spectral gain at diurnal frequency.

Section 5 develops, in contrast, a somewhat more sophisticated implementation over "short time scales" with a strong non stationary trend in the pore pressure signal. This occurs in the case of a moving excavation front passing nearby the pressure sensor. The objective is to identify the effect of the evolving "EdZ" at the location of the pressure sensor, and in particular, its effect on the evolution of Ss(t). The time scale of analysis is typically a couple of weeks, and a new tool has been developed (based on statistical envelopes) in order to identify the modulation in time of the amplitude of earth tide fluctuations.

Section 6 provides a summary of the methods and results, and suggests possible modifications and extensions of this work concerning statistical methods (multi-cross analyses) and hydro-geologic interpretations (enhanced hydro-mechanical models and inverse problems): see the "outlook" in Sub-Section 6.2.

The text is completed by an appendix (for abbreviations and symbols), and by a list of references.

Hydrogeologic setting and pressure signal measurements

Hydrogeologic setting at Mont Terri (Opalinus clay rock)

The Mont Terri Rock Laboratory was excavated within the "Opalinus Clay" of the Mont-Terri anticline, which has a very low permeability. The geological cross section shown in

Figure 1 shows the asymmetrical, dome-shaped fold of the Mont Terri anticline, where sediments of Triassic and Jurassic age are exposed.

The Opalinus clay is about 180 million years old (Mesozoic age, Jurassic period, Lower Dogger, Aalenian). It is a shallow marine mud deposit, extending over wide areas in Central Europe. It is underlain by the Jurensis Marls (Toarcian) and overlain by the "Lower Dogger" limestones (Bajocian). Three facies have been identified within this clay rock: sandy, carbonate-rich sandy, and shaly. The latter facies is studied as an analogue for hypothetical disposal of radioactive waste in a naturally isolating geologic formation (the Mont Terri site is a research site only, not destined to host nuclear wastes or any radioactive products whatsoever). For more details on geology and hydrogeology of the site: see [START_REF] Fatmi | Statistical pre-processing and analyses of hydrogeo-meteorological time series in a geologic clay site (methodology and first results for Mont Terri's PP experiment)[END_REF], [START_REF] Fatmi | Méthodologie d'analyse des signaux et caractérisation hydrogéologique: Application aux chroniques de données obtenues aux laboratoires souterrains du Mont Terri[END_REF] and [START_REF] Schaeren | Tunnels du Mont Terri et du Mont Russelin. La traversée des "roches à risques": marnes et marnes à anhydrite[END_REF]. See also, more recently, [START_REF] Nussbaum | Analysis of tectonic structures and excavation induced fractures in the Opalinus Clay, Mont Terri underground rock laboratory (Switzerland)[END_REF].

We emphasize again that the Opalinus clay rock at Mont Terri has very low permeability or hydraulic conductivity: K ≈ 10 -13 to 5 10 -13 m/s (1 E-13 to 5 E-13 m/s), ± 1 order of magnitude, in the intact zone outside the Excavated Damaged Zone. This clay rock also displays a significant self-sealing capacity due to a high proportion of swelling minerals (about 10%). In contact with water, it swells and tends to seal the cracks and fractures that occurred after excavations (galleries, niches). This self-sealing reduces permeability and participates in restoring initial rock properties.

These factors, among others, explain the choice of this type of clay rock formation by the Swiss authorities -and by the international Mont Terri consortium -as a potential host rock for high level radioactive wastes. The Mont Terri URL (URL stands for Underground Research Laboratory ) is managed by SWISSTOPO. The Mont Terri consortium includes Swiss, French, British, Canadian, Spanish, Belgian, U.S., Japanese, and other nations' public authorities and their national research institutions.

Other types of geological formations were also studied by the Swiss authorities, including a crystalline rock (granite) at the Grimsel site in the Swiss Alps (e.g., FEBEX drift in situ experiment). Thermo-Hydro-Mechanical upscaling and modelling of the FEBEX drift in situ experiment was studied in Cañamón's thesis (2006), later published as [START_REF] Valera | Coupled phenomena in 3D fractured media (analysis and modeling of thermo-hydro-mechanical processes[END_REF]. These references also include statistical signal analyses of the FEBEX thermo-migration laboratory experiment, or "mock-up" experiment (see also [START_REF] Cañamón | Wavelets and statistical techniques for data analysis in a mock-up high-level waste storage experiment[END_REF]. The in situ and mockup FEBEX experiments were performed under funding by the Spanish nuclear authorities ENRESA.

In the remainder of this text, we focus exclusively on the analysis of pressure signals collected in situ in Mont Terri's Opalinus Clay rock site.

Pressure signal measurements and pre-processing (Mont Terri BPP-1 borehole)

The pressure signals used in this work were obtained over a period of 10 years with halfhourly time steps (on average), at the Mont Terri URL, in the framework of the LP14 project (Long term Pressure monitoring project -phase 14). Concerning Mont Terri's LP project, see [START_REF] Marschall | Chap. 8: Hydrogeology experiments[END_REF]. The BPP-1 borehole (101 mm diameter, and 20 m length) was dug in the year 1996, from niche PP located in Mont Terri's "exploration gallery". The measurement sections or sensors (PP-1 and PP-2) are located at the end of the BPP-1 borehole, outside the excavation damaged zone (EDZ) of the exploration gallery, in the siltyshaly facies of the Opalinus clay. Absolute pore pressures P(t) were measured from boreholes (as explained further below), and atmospheric pressure P AIR (t) or P ATM (t) was measured underground in the exploration gallery. Absolute pore pressures were then transformed into relative pressures, defined as: P REL (t) = P(t)-P ATM (t). For more technical details on the BPP-1 borehole and its pressure sensors, see Figure 3 further below, and see [START_REF] Thury | Mont Terri rock laboratory. Results of hydrogeological, Geochemical and geothechnical Experiments Performed in 1996 and 1997[END_REF]. Finally, see also the previous Figure 1 for 3D perspective views of the tunnels, galleries, and boreholes.

Signal pre-processing, record lengths, missing data and other issues

The BPP-1 borehole was initially selected for this study because it provided the longest pore pressure times series over a period of about 10 years (from 17/12/1996 to 30/06/2005). However, after pre-processing for outliers, data gaps, and variable steps, the 10 year long sequences from the PP1 and PP2 pore pressure sensors (and from the air pressure sensor) had to be fractioned into shorter "reconstituted" contiguous sequences. For relative pressure P REL (t), the resulting sequences were roughly on the order of a year (no longer than 1.5 year at best).

It should be emphasized that, in order to obtain relative pressure sequences (P REL (t) = P(t)-P ATM (t)), the pore pressure P(t) and air pressure P ATM (t) had to be processed "jointly" [START_REF] Fatmi | Statistical pre-processing and analyses of hydrogeo-meteorological time series in a geologic clay site (methodology and first results for Mont Terri's PP experiment)[END_REF] Statistical tests, involving recursions, were used to detect possible outliers (transformed into missing data). In the presence of missing data, for partial reconstitutions, we used:

decompositions into trends and residuals using moving average filters (and other filters), -linear interpolation techniques (usually only for relatively short data gaps), -a bidirectional first order autoregressive model (going forward and backward in time) -statistical tests and validations of reconstituted series based on moments and covariances Time step homogenization is also a kind of partial reconstitution. We used for this a linear interpolation/extrapolation method, with a constant time step "Δto" selected by the user.

Here, all processed signals have, after "homogenization", a half-hourly time step (Δto = 30 mn). This choice is in agreement with the nominal acquisition step (also the most frequently encountered time step in the collected data).

Figure 2 shows a simplified flow chart presenting these Pre-Processing tasks and methods.

For more details on pre-processing methods, see [START_REF] Fatmi | Méthodologie d'analyse des signaux et caractérisation hydrogéologique: Application aux chroniques de données obtenues aux laboratoires souterrains du Mont Terri[END_REF], Chap.3) and Fatmi et al. (2007[START_REF] Ababou | Quantitative Stochastic Hydrogeology: the Heterogeneous Environment. Chap. 8 in Part III of "Overexploitation & Contamination of Shared Groundwater Resources: Management, (Bio)Technological, and Political Approaches to Avoid Conflicts[END_REF]. At the end of the pre-processing stage, the processed (transformed) signal is ready for statistical analyses. The processed signal may be shorter than the initial "raw" signal.

Pressure sensors location with respect to the excavated Gallery 1998

Since we are especially interested here in the effects of excavation on the hydrogeologic properties of the clay rock (EDZ), let us focus here on a brief description of the excavation of gallery Ga98.

The Ga98 gallery was dug at Mont Terri during the year 1998, over a period of 5 months (more precisely, from 17/11/1997 to 23/04/1998).

Recall that the main goal of this work is to study not only the properties of the claystone, but also, their evolution in response to excavation (Excavation Damaged Zone). This will be done, for the case of gallery Ga98, by analysing the effects of excavation on pore pressures collected at sensors PP1 and PP2, knowing that they are distant only a few meters (4.4 m and 4.5 m) from the excavated gallery Ga98 (Figure 3).

Thus, due to their close location to the 1998 gallery, the sensors of borehole BPP-1 can be used to monitor pore pressures during the excavation of the 1998 gallery: we will exploit this particular configuration in this paper, as in a "mine-by test".

Note on related analyses of other signals at Mont Terri and at other clay sites

While we focus in this work on air and pore pressures at Mont Terri, for the sake of completeness, we present in this section a brief note on other signal analyses: e.g., air temperature, relative humidity, and crack aperture at Mont Terri, and also, signals recorded at other URL sites in clay rock (Bure and Tournemire, in addition to Mont Terri). Thus:

The Mont Terri "BPP-1" pore pressure signals were compared and cross-analysed with air pressure, air temperature, and relative air humidity measured in the exploration gallery of the Mont Terri URL.

In addition, similar signals were analysed and cross-analysed at two other clay rock URL's (Underground Research Laboratories) in addition to the Mont Terri URL:

• the clay rock URL of Tournemire (Aveyron, France), operated by the French "Institut de Radioprotection et de Sûreté Nucléaire" (IRSN), and • the clay rock URL of Bure (Meuse / Haute Marne, France), operated by the French national agency on the management of nuclear waste (ANDRA).

Concerning these sites (Tournemire and Bure), the reader is referred to the Ph.D. thesis of [START_REF] Fatmi | Méthodologie d'analyse des signaux et caractérisation hydrogéologique: Application aux chroniques de données obtenues aux laboratoires souterrains du Mont Terri[END_REF], available on line. Also, it should be noted that despite their similarities, the three clay rock sites have different tectonic, structural and hydro-mechanical characteristics (e.g., see Nussbaum et al. 2007, comparing Bure and Mont Terri).

• Furthermore, crack aperture time series measurements were also developed at the Mont Terri site (Cyclic Deformation phase 14 project "CD14"). The crack displacement signal was statistically analysed. The air temperature, pressure, and relative humidity were also measured and cross-analysed together with crack displacements. Some of the statistical methods used in the present study (cf. Section 3) were also applied to the "cross-analysis" of crack aperture response to hydro-meteorological fluctuations at various time scales. The results of these crack aperture investigations are reported in Möri et al. (2012, in press at the time of this writing).

Gallery 1998

Fig. 3. Tunnels, galleries, and boreholes at the Mont Terri URL, showing Gallery 1998 ("Ga98"), and borehole BPP-1 dug from the PP niche. This niche is located in the reconnaissance gallery ("exploration gallery") parallel to the main road tunnel (dark grey).

The pore pressure sensors used in this study are located near the end of borehole BPP-1 (lower left corner of blue rectangular frame), only 4 m away from Ga98. Their response during excavation of Ga98 is studied in this work.

Previous results and scope of statistical analyses this work (long and short time scales analyses of pore pressures)

The general objective of this work, in the context of radioactive waste disposal in geologic clay sites, was to exploit pressure measurements in order to quantify the phenomena affecting pore pressure (earth tides, barometric pressures,…), and in fine, to estimate the hydro-geologic or hydro-mechanical properties of the clay rock : elastic compressibility and/or specific storativity (Ss), effective elastic porosity(Φ), ….

Some pore pressure analyses methods and results (Fatmi et al. 2007[START_REF] Fatmi | Statistical pre-processing and analyses of hydrogeo-meteorological time series in a geologic clay site (methodology and first results for Mont Terri's PP experiment)[END_REF] All these analyses, as mentioned above, exploit the effects of natural perturbations of pore pressure, caused by hydro-meteorological as well as geophysical oscillations:

• hydrometeorology: atmospheric pressure oscillations; • planetary geophysics: earth tides.

Thus, in the short time scale analyses (excavation works) we consider also the effects of changes in hydro-geological properties due to human activities in the URL site. The excavation of an underground gallery serves as a "mine by" test. The objective is to monitor the response of pore pressure (measured at a fixed location) to the evolution of the excavation front. The effect of the excavation (if any) should be superimposed on the natural perturbations of pore pressure (atmospheric and earth tide effects).

We will see that it is possible to detect, statistically, the effects of the evolving excavation front on the component of pore pressure fluctuations due to earth tides. The time scale of analysis here is relatively short, because we focus on a particular "non stationary" event, the approach of the excavation front close enough to the pressure sensor(s). This is to be distinguished from the "long time sale analyses", where we try to use records that are as long as possible, and not too much influenced by local trends.

Along these lines, for short time scale "excavation effects", our main objective is as follows:

1. to analyse the effects of drift excavation on pore pressures during excavation (synexcavation)1 , and if possible, to analyse this effect continuously as a function of time (statistically, this boils down to analysing a signal with strongly non stationary fluctuations).

In the case at hand, we focus on the passage of an excavation front of gallery GA98 at Mont Terri (gallery 1998), and we examine closely its possible influence on the pore pressures measured in pressure sensors (PP1, PP2) of the BPP-1 borehole (here, we will narrow this down to the single pressure sensor PP1, closest to the excavation front). The specific objectives are, essentially:

• to quantify the evolution of elastic specific storativity Ss(t) as a function of time, during the different phases (pre-, syn-, post-excavation phases), and partiuclarly, the maximum change of Ss(t) during the syn-excavation phase; • to interpret the evolution of Ss(t) and to compare it with the distance D(t) of the moving excavation front (and the moving "EdZ") relative to the location of the fixed pressure sensor.

Methods for statistical analyses and interpretations of pressure signals

Pre-processing procedures were described in the previous introductory sections (see flowchart in Figure 2). We now focus on methods of statistical analyses, and also, on the simplified hydro-mechanical "models" to be used for hydro-geologic interpretation of pore & air pressure analyses.

Overview of methods of statistical analysis (and interpretation)

After the preliminary stage of signal pre-processing, we are left with several records of contiguous time series with constant time steps (some of these signals have been preprocessed to some extent, and possibly, they have been partially reconstituted). These preprocessed pressure time series are then statistically analysed in order to characterize the hydraulic and hydro-mechanical behaviour of the porous clay rock formation (Ss, Φ). To achieve this, the following signal analysis methods were used and applied to Mont Terri signals2 :

i. Auto-correlation and cross-correlation analyses of pressure signals versus time lag, and also (closely related), deconvolution or identification of input/output transfer function, e.g., between atmospheric pressure and relative pore pressure. ii. Fourier spectral analyses in the frequency domain, based on Fourier transforms of the correlation functions. This leads to identification of dominant frequencies, and also, cross-spectrum analysis and identification of the spectral gain between two signals (such as atmospheric pressure and relative pressure). iii. Orthogonal multi-resolution wavelet analysis was used to decompose the signals in the time/scale domain (two-parameter domain). Thus, for any selected dyadic time scale (e.g., for the semi-diurnal time scale of earth tides), we isolated and extracted the corresponding "wavelet component" of the original signal. Also, we used wavelets in some cases as a filter, decomposing the original signal into a filtered quantity (the wavelet "approximation") and its residual signal (the wavelet residual) at any given dyadic time scale (the scale of the filter). iv. Estimation of the statistical envelope of a modulated random signal based on Cramer-Leadbetter statistical envelope theory, implemented with the Hilbert transform. This technique is used to characterize the modulation, in time, of the signal extracted from semi-diurnal wavelet component analysis of pore pressure during excavation.

Each of the above-mentioned techniques was used at one stage or another in our study of Mont Terri pressure signals on various time scales.

The correlation and cross-correlation function analyses versus time lag (I), as well as the corresponding Fourier spectral analyses in frequency space (II), are both based on the theory of random processes: the reader is referred to [START_REF] Bras | Random Functions in Hydrology[END_REF]; [START_REF] Box | Time Series Analysis, Forecasting, and Control[END_REF]; [START_REF] Papoulis | Probability, Random Variables, and Stochastic Processes[END_REF], [START_REF] Yaglom | Correlation Theory of Stationary and Related Random Functions: Basic Results[END_REF], [START_REF] Yevjevich | Stochastic Processes in Hydrology[END_REF]. Concerning spectral estimation and analysis (II), see also: [START_REF] Blackman | The Measurement of Power Spectra[END_REF], [START_REF] Max | Méthodes et techniques de traitement du signal et application aux mesures physiques[END_REF], [START_REF] Lacoume | Méthodes et techniques de traitement du signal[END_REF], [START_REF] Priestley | Spectral analysis and time series[END_REF]), Vanmarcke (1983). Time domain deconvolution and input/output analyses, causal or non-causal, linear or non linear, are presented in Labat et al. (1999c;2000a) for the case of rainfall-runoff signals.

Multi-resolution wavelet analysis (III) is an orthogonal representation of signals in "scale-time" space (a two-parameter space). The mathematical theory is presented in [START_REF] Mallat | A Theory For Multiresolution Signal Decomposition: The Wavelet Representation[END_REF]. See Labat, Ababou et al. (1999a,b;2000b;2002) concerning applications of wavelets to rainfall/runoff signals in karst hydrology. In multi-resolution wavelet analyses, "time" is discrete in the usual way (t(n) = n×Δt), but wavelet "scales" are dyadic (powers of 2 : 2 m ×Δt).

Wavelet "scales" are analogous to inverse frequencies, but still, they should not be confused with Fourier frequencies.

The theoretical tools III (wavelets) and IV (envelopes) will be combined and used together in Section 5 (5.2+5.3). The goal is to assess, on a short time scale of a few weeks, the influence of excavation on the evolution of pore pressure and, finally, of clay rock properties in the "EdZ".

The mathematical and statistical methods cited above have well known theoretical bases, however, they can be implemented in different ways in practice. For this study, we have devised an integrated set of programs, written in MATLAB and organized in "TOOLBOXES" as follows:

• Signal Pre-Processing TOOLBOX, and • Signal Analyses TOOLBOX.

The two TOOLBOXES are described in Fatmi 2009's PhD thesis, available on line.

3

It should be noted that some of the methods and techniques were specifically designed or adapted by us for this study, e.g.: estimations of cross-spectra and spectral gains; identification of time domain transfer functions; filtering techniques (comparisons between moving averages and wavelet approximations); and also, an original implementation of the Cramer-Leadbetter envelope, performed with the Hilbert Transform, and applied to a selected wavelet component.

Other auxiliary studies not reported here include the design of a number of synthetic tests to validate statistical methods (cross-correlations, cross-spectra, envelopes, comparison of filtering techniques, etc). For details on some of these tests, the reader is referred to [START_REF] Fatmi | Méthodologie d'analyse des signaux et caractérisation hydrogéologique: Application aux chroniques de données obtenues aux laboratoires souterrains du Mont Terri[END_REF].

Figure 4 shows a simplified flowchart presenting some of the signal analysis methods used in this work. These methods were programmed in a MATLAB TOOLBOX for SIGNAL ANALYSES (distinct from the PRE-PROCESSING TOOLBOX described in Figure 2).

Univariate (auto) correlation functions, bivariate (cross) correlation functions, and temporal transfer functions

Correlation function analysis is a way to study the temporal structure of the signals, in the time lag domain. The univariate analysis of a single signal leads to a study of its autocorrelation function. The bivariate analysis of two signals jointly leads to both auto-and cross-correlation functions. More precisely, the temporal structure of a signal X(t), or of two signals X(t) and Y(t) jointly, "is captured" by a temporal structure function that depends on two points in time (t, t+τ).

The stationarity hypothesis simplifies this analysis by assuming that the structure function depends only on the time lag (τ). In addition, the ergodicity hypothesis (which is notoriously hard to validate in practice) further assumes that the time average of the signal is equivalent to the ensemble mean (mathematical expectation) as the record length becomes infinitely long.

Thus, for two statistically stationary and ergodic processes, X(t) and Y(t), the crosscovariance function is supposed to depend only on the time lag (τ) between X(t) and Y(t+τ) (and it does not depend on "t"). The cross-covariance Cxy (and the cross-correlation Rxy) can be estimated as follows, from a single data record of the two processes (X(t),Y(t))

4 : 1 1 ( ) ( ) ( ) i N j XY i i j i C j X t X Y t Y N = - - - + =     = - -          and ( ) ( ) . XY XY X Y C j R j σ σ = (1)
In this equation:

• "i" is the discrete time (t(i)) • "j" is the discrete time lag (τ(j)) • R XY (j) is the cross-correlation (function of discrete time lag "j"), • C XY (j)
is the cross-covariance (function of discrete time lag "j"), •

N: total number of data in the observed time series.

One also defines "M", the maximum time lag to be retained for correlation function analyses.

The discrete time lag index "j" runs through j = ±{0,1,2,…, M}. It is necessary to keep M ≤ N-1

(in practice M ≤ N/3 or much less).

A sub-sampling step k 0 can be defined, with k 0 = 1 or 2 or 3 or ... M (take k 0 = 1 if "no subsampling"). The time lag in physical time units is: τ j = j.Δt (if k 0 = 1) or more generally τ j = j.k 0 .Δt (k 0 ≥ 1) where ko is the sub-sampling step. In our case, we choose the obvious, i.e., the smallest sub-sampling step (k 0 =1) in order to sweep all the samples of the time series. For very long series or unduly small time steps, sub-sampling (k 0 ≥ 2) may be advisable, but this is not our case here.

The above Eq. 1 gives a biased estimate of the cross-covariance function Cxy(τ). To obtain an unbiased estimate, the denominator (N) should be replaced by (N-j-1) in the above estimator, although some authors also use (N-j).

Finally, it is noted that Cxy(τ) is directly related to (but not equal to) the temporal transfer function Fxy(τ) between processes X(t) and Y(t). See Labat et al. (1999c;2000a) for a detailed analysis and discussion of temporal transfer functions in hydrogeology & hydrology (Instantaneous Unit Hydrographs). Also, it is interesting to note that the transfer function Fxy(τ) can be interpreted like the "slope" of a generalized linear regression between the two signals X(t) and Y(t). The identification of this transfer function in time is an inverse problem; it is also a particular type of optimal estimation problem, also known as the "deconvolution problem".

Spectral and cross-spectral analyses, frequency gain

Fourier spectral analysis decomposes the signal (or its structure function) into periodic functions, expressed in the frequency domain using the Fourier Transform: e.g. the Fast Fourier Transform (FFT) implemented in MATLAB and many other software's.

More precisely, univariate spectral analysis consists in studying the frequency spectrum of only one signal (its auto-spectrum S XX (f)), while bivariate or cross-spectral analysis, treats simultaneously two signals, leading to the frequency cross-spectrum S XY (f) of the two signals.

Other related quantities are the co-spectrum and the phase spectrum (not used here), and the frequency gain (used below). The frequency gain can be interpreted as a spectral transfer function, particularly in the case where one of the two signals (X(t)) is the "cause" of the other signal (Y(t)). This frequency gain is (in a sense) the spectral version of the temporal transfer function between the "input" X(t) and the "output" Y(t).

The spectrum Sxx(f) of a signal X(t) can be defined and estimated in various ways (continuous/discrete time, infinite/finite time window, with/without filtering). Below, the spectrum of X(t) is expressed as a function of frequency, for the case of a discrete time process, within a limited time window, and with filtering…

Estimation of the Fourier spectrum (Wiener-Kinchine method, filtered)

The Fourier auto-spectrum of signal X(t) is estimated here as the Fourier transform of the correlation function Rxx(τ), based on the Wiener-Khinchine theorem, valid for stationary random processes. The final expression proposed here includes, also, a time-lag filter "D(τ)".

The dimensional estimated spectrum S(f) is then of the form:

( )

1 ( ) 2 (0) 2 . ( ).cos 2 M XX i XX j XX i j S f t C D C j j t f π =   = Δ + Δ        (2a)
and the reduced dimensionless spectrum s(f) is of the form:

( )

1 ˆ( ) 2 1 2 . ( ).cos 2 M XX i j XX i j s f D R j j f π =   = +        (2b)
where:

2 ( ) ( ) XX XX XX R j C j σ = , ( (0) 1)
XX R = . (3a)
C XX (τ) is the auto-covariance function and R XX (τ) is the auto-correlation function, or normalized covariance function. Similarly, S XX (f) is the auto-spectrum, and "s XX " is the reduced dimensionless spectrum (the reduced spectrum is normalized by the variance of the signal and expressed in terms of dimensionless frequency i f ˆ). Furthermore, in order to obtain a frequency dependent "measure" of fluctuations analogous to the notion of standard deviation, we define the "Root Mean Square spectrum" as follows (it has the same units as the signal X(t) and σ X ):

( ) 2 ( )

RMS XX i XX i S f S f f = Δ . ( 3b
)
The frequencies and dimensionless frequencies are defined as follows:

i f : frequency (Herz) (1 Hz=1 s-1); 1 i i T f = : period (seconds)
Let M be either the number of time steps in the data, or the maximum number of lags allowed in the Tuckey filter D(j) (see next section 3.3.2). If the Tuckey filter is turned off, one can let M = N-1 where N is the total number of points in the observed signal. Then, we have:

1 2 i i f M t = Δ : period (seconds)
According to Shannon's sampling theorem (see Figure 5), here are the maximum and minimum frequencies that can be "identified" given the time step and the duration of the signal5 :

max min 1 1 2 f Δt T = = [Hz] (4a) min max max 1 1 1 2 2 f T t M t = = = Δ [Hz], (4b)
where max t M t = Δ (size of the time window). In these relations, Tmin and Tmax are, respectively, the smallest and largest periods that can be represented with this time series

( max min 2 T T ≥
). The dimensionless frequency ˆi f is defined by:

ˆi i f f t = ×Δ  0 ˆ2 i i f M k = × × ; ( 4c
)
where ko is the sub-sampling step in time (we always use full sampling with ko = 1 in this work). The dimensionless frequency ˆi f belongs to the interval [0, 1/2], that is:

0 0 . 5 i f ≤ ≤
. Thus, for k 0 = 1, we have:

0 ˆ0 f = and ˆ0.5 M f = (4d)
The index "j" in the Fourier spectrum/covariance Eq. 2 is the label of the time lag, defined as follows:

j = 0, 1, 2, …, M ≤ N-1
where the zero time lag (j = 0) is taken into account separately outside the sum in Eq. 2. On the other hand, the index "i" is the label of discrete frequency, or of dimensionless frequency

( ˆi f ): i = 0, 1, …, M ≤ N-1 or (see the note just below): i = 0, 1, …, M, …, N-1.
The number "M" represents the discrete size of the filter truncation window (number of Δt's in the maximum lag considered), while "N" is the total number of data of the observed signal (full time window), and Δt is the constant time step of measurements (possibly after homogenization of time steps). In general, for a filtered spectrum, M ≪ N (at most, M=N-1 for the unfiltered spectrum).

Note. Concerning the discrete number of frequencies represented in the filtered spectrum (index "i"), it is possible to impose either a "coarsened" representation of the spectrum with frequencies running from i=1 to i=M (where usually M≪N), or alternatively, to maintain a fine representation of the spectrum with frequencies running from i = 0, 1, …N-1, regardless of the maximum time lag "M" used in the Tuckey filter "D". In the latter case, one should replace "M" by "N-1" in all of Eqs. 4 (index "i").

The Tuckey filter ("hanning" filter) for spectral estimation

Finally, note that "s(f)" in Eq. 2 is an estimator of the spectrum, with a filter function "D". More precisely, "Dj" represents a time-lag filter function; it is applied to the correlation function within the Fourier transform. The weight function D(j) included inside the Fourier transform represents a "smoothing function" or "filter" of Rxx(j) in time lag space (j), where "j" is the discrete index of the dimensional lag time τ(j). Thus, we can write equivalently D(j) or D(τ(j)).

Our particular choice here the "Tukey-2" filter (also called a "hanning" filter):

1 cos 2 j j M D π     +         = (5)
where the index "j" runs generally over (j = -M,...,-1,0,+1,...,+M), or here over (j = +1,...,+M), in the particular discrete Fourier Transform sum written in Eq. 2.

The filter D(j) is necessary in most cases. Indeed, it is necessary to apply smoothing to Rxx(τ) or Rxx(j) because of the sampling errors on the estimated Rxx(j), which is obtained on a finite time window. Thus, it can be shown precisely that the estimation error on the covariance Cxx(j) increases with time lag τ(j): see Yevjevich (1972, Chap.3) and Bartlett (1946). For example, in the case of an autoregressive process with exponential covariance function

( ) ( ) 2 exp XX X COR C τ σ τ τ = -
, it has been shown (Ababou et al. 1988: Appendix 2) that the variance of sampling error, defined by

( ) ( ) ( ) ( ) 2 2 ˆXX RMS XX Var C C ε τ τ τ =
, is given by:

( ) . This illustrates why it is necessary to apply a filter like D(τ(j)) in order to smooth out the estimated Cxx(τ(j)) or Rxx(τ(j)) at large time lags, and also to smooth out the resulting frequency spectra obtained by Fourier transform of Cxx(τ(j)) or Rxx(τ(j)).

2 1 1 1 2 1 1 2 4 2 COR COR RMS MAX MAX MAX COR MAX t t t t τ τ τ τ τ ε τ τ     = - - + -           . ( 6
The filter we have chosen here ("Tukey-2", or "hanning" filter) adequately smoothes out large time scales (or noise at the scale of low frequency increments in the frequency domain). Experience with hydrological signals indicates that this filter has less variance bias than others (according to [START_REF] Mangin | Pour une meilleure connaissance des systèmes hydrologiques à partir des analyses corrélatoires et spectrales[END_REF], it over-estimates the total variance by no more than 8%). However, other filters may be considered in future: in particular, it should be noted that a "minimum bias" spectral filter has been proposed by [START_REF] Papoulis | Probability, Random Variables, and Stochastic Processes[END_REF].

Finally, note that the Fourier transform of our time-lag filter D(τ(j)) is itself a "filter" in frequency space (D(f)). Its spectral band should contain a large part of the spectral band of the signal, in order not to waste too much information in this filtering procedure when estimating the spectrum. In our applications, the maximum lag "M" in the D(τ(j)) filter (Eq. 5) was adjusted depending on each specific application and objective (time scales and frequencies of interest).

Cross-spectrum, frequency gain, and coherence function

The (dimensional) cross-spectrum of two signals (X(t), Y(t)) is a complex function of frequency which can be described as follows:

( ) ( ) ( ) XY i f XY i XY i S f S f e ϕ = or Re Im ( ) ( ) ( ) XY i i i S f S f S f = + (7a)
where

2 2 Re Im ( ) ( ) ( ) XY i i i S f S f S f = + (7b) and 1 Im Re ( ) ( ) tan ( ) i XY i i S f f S f ϕ -  =     , (7c)
are respectively the amplitude (cross)-spectrum and the phase (cross)-spectrum. The latter is defined only in an interval of size π, such as [-π/2,+π/2]. More precisely, however, the phase spectrum can be identified unambiguously within the full interval [-π, +π] with the equation:

( ) 1 Im ( ) ( ) ( ) cos ( ) R i XY i i XY i S f f sign S f S f ϕ -  =     , (7d)
as will be demonstrated shortly with a simple test.

Given the complex cross-spectrum Sxy(f), one can also define its real part (or "co-spectrum"), and its imaginary part (or "quadrature spectrum"). These are defined and estimated using the Wiener-Khinchine representation (similar to that used in Sec. 3.3.1 for the autospectrum). The result is expressed below in terms of sine and cosine transforms (where a Tuckey filter "Dj" is inserted):

( ) Re 1 ( ) 2 (0) ( ) ( ) cos(2 ) M i XY X Y X Y Y X j i j S f t R R j R j D j t f σ σ π =   = Δ + + Δ        (8a) ( ) Im 1 ( ) 2 ( ) ( ) sin(2 ) M i X Y X Y Y X j i j S f t R j R j D j t f σ σ π =   = Δ - Δ        (8b)
It is useful to define also the dimensional Root Mean Square cross amplitude spectrum, as follows (it is important to note that this RMS cross-spectrum has the units of σ X ×σ Y ):

( ) 2 ( )

RMS XY i XY i S f S f f = Δ (9)
Finally, we define a spectral coherence function:

( ) ( ) ( ) ( ) RMS XY i RMS XY i RMS RMS XX i YY i S f Coh f S f S f = (10)
and a spectral gain, or frequency gain function:

( ) ( ) ( ) RMS XY i RMS XY i RMS XX i S f G f S f = (11)
In the latter equation defining the "gain", it should be noted that the two signals are treated differently: X(t) is considered the input signal, and Y(t) the output signal. Thus, Gxy(f) is the input→output spectral gain. It quantifies, for each frequency "f", the effect G of the signal X(t) on the output signal Y(t). This effect is called the gain.

Synthetic example of cross-spectral analysis on two harmonic signals

Two "synthetic" harmonic signals are used in order to test the auto-and cross-spectral estimates described above (spectrum of each signal, then cross amplitude spectrum, phase spectrum, frequency gain). 6 The two signals, input X(t)=X 1 (t) and output Y(t)=X 2 (t), are defined as follow:

4 1 0 1 ( ) cos( ) k k k k X t A A t ω φ = = + +  ; 4 2 0 1 ( ) cos( ) k k k k X t B B t ω γ = = + +  (12a)
Both signals have four harmonic components, with the same four frequencies ω i = 2π/T i (i = 1,2,3,4) with periods T i = 6 months, 8 days, 24 hours and 12 hours, respectively. In addition, each signal has its own constant mean (A 0 or B 0 ). The mean may be viewed as the zero frequency harmonic (i = 0). The four amplitudes and phases are different for the two signals.

Here is a summary or the values used in this test:

-Signal X 1 (t):

Mean A 0 = 50 [units of X 1 ]; amplitudes (A i ) = (12, 9, 4, 1) [units of X 1 ]; phases (φ i ) = (π/4, 0, π/8, π/6) [rad]. - Signal X 2 (t): Mean B 0 = 50 [units of X 2 ]; amplitudes (B i ) = (1, 4, 9, 12) [units of X 2 ]; phases (γ) = (0, 0,0, 0) [rad].
In the test presented below, the two synthetic signals are generated during a time span of 1 year (in comparison, the longest period in the signal is 6 months). The time step is Δt = 900 s = 15 mn (in comparison, the shortest period in the signal is 12 hours). Spectral estimations are carried out "unfiltered" (the Tukey filter "Dj" is not applied: we take M = N-1 as the maximum time lag).

The spectral estimation results, shown in Figure 5 & Figure 6, are exactly as expected: they coïncide very closely with the exact theoretical values that we have inferred directly from cross-spectral analysis of these harmonic signals (assuming infinite record t MAX →∞, and continuous time Δt→0). For reference, here are some theoretical results for each discrete frequency f(i)=ω i /2π:

( ) 1; ( ) ; ( ) / ; ( ) RMS RMS XY k XY k k k RMS XY k k k XY k k k Coh f S f A B G f B A f ϕ φ γ = = = = - (12b)
For instance, focusing on the frequency of the 24 hour harmonic, we have the following results (we use unbiased covariances, except for phase & gain spectra which are estimated using biased covariances):

• the RMS auto-spectrum of X 1 (t) should be exactly 4.00 [units of X 1 ], and we find 3.996 from the estimated spectrum (Fig. 7.a); • the RMS auto-spectrum of X 2 (t) should be exactly 9.00 [units of X 2 ], and we find 8.999 for the estimated spectrum (Fig. 7.b);

• the RMS coherence spectrum [dimensionless] should be exactly 1 for each of the four harmonics, and we do find that result for the estimated coherence spectrum (Fig. 7.c); • the RMS amplitude cross-spectrum should be exactly 6.00 [units of √X 1 √X 2 ], and we find 5.997 for the estimated cross spectrum amplitude (Fig. 7.d); • the phase spectrum should be exactly π/8 [rad], and we find 0.391 [rad] for the estimated phase (Fig. 6.a); • the spectral gain between X 1 (t) and X 2 (t) should be exactly √9/4 = 3/2 [units of √X 2 /√X 1 ], and we find 1.50, the exact result with two significant digits (Fig. 6.b).

As a final note on spectral estimation, it should be emphasized that the Wiener-Kinchine method used in this work relies on the Fourier transform of auto-and cross-correlation functions, which are themselves estimated in the time domain assuming that the processes (X(t), Y(t)) are jointly statistically stationary. Alternatively, however, the power spectrum Sxx(ω) of a signal X(t) could be estimated from the "periodogram" Pxx(ω), which is the direct Fourier transform of the signal X(t i ) itself (over a finite time window, in discrete time).

Filtering or smoothing will be again necessary, as in the Wiener-Khinchine method. The difference is that the required smoothing can now be performed in frequency space (ω). For instance, the periodogram can be smoothed out using a moving average filter in frequency space, to finally obtain a filtered estimate of Sxx(ω). For more details on this, the reader may consult Blackman and Tuckey (1958), Yevjevich (1972, Chap.3), [START_REF] Papoulis | Probability, Random Variables, and Stochastic Processes[END_REF]Pillai (2002, Chap.12), and[START_REF] Priestley | Spectral analysis and time series[END_REF].

Multiresolution wavelet analysis

Introduction on wavelets

The analysis of a signal in terms of multi-resolution wavelets uses an orthogonal basis

( ) { } , j k t ψ
where "k" represents the temporal translation (position), and "j" represents the dilation or the contraction of time scale, based on a dyadic time discretization of the signal X(t).

The basis being orthogonal, the components of different "scales" are independent from each other. Moreover, each component at a given scale can be traced according to the "time" parameter or "temporal position". This is to be compared with the traditional Fourier analysis, where each Fourier coefficient C(f) represents the overall contents of the signal at frequency « f ». In other words, Fourier coefficients C(f) encapsulate the "f" frequency content of the signal X(t) for all times "t", without representing the possible temporal localization of the frequency content.

On the other hand, the orthogonal wavelet approach makes it possible to clearly highlight non-stationary behaviour and time-localization phenomena (Labat et al. 2000b). It also makes it possible to follow the relationship between the various scale levels as a function of time, and to isolate each independent component to study its own evolution. In this sense, these analyses can lead to a better identification of the process responsible for the observed variations. For an extension to wavelet multivariate analyses, or cross-wavelet analyses, see [START_REF] Labat | Analyse multirésolution croisée de pluies et débits de sources karstiques[END_REF].

We choose here to conduct multi-resolution analysis with discrete orthogonal wavelets, using the wavelet bases of Daubechies. The software employed in this work combines a set of wavelet functions programmed in MATLAB and available on line at the University of Stanford (WAVELAB), and other functions contained in MATLAB's own WAVELET TOOLBOX (some of which are used by Stanford's WAVELAB package).

Wavelet decomposition of a signal on an orthogonal multiresolution basis

In this subsection, we present the decomposition of a signal on a wavelet basis, first for continuous wavelets (non orthogonal bases), and secondly for dyadic multi-resolution wavelets. In the remainder of this text, we will focus essentially on the latter wavelets, for which a compactly supported orthogonal wavelet basis can be constructed (Daubechies basis). The term "dyadic" refers to the fact the scales of multiresolution wavelets vary discretely as powers of "2". For example, if the length of the recorded signal is t MAX = 64Δt, then there are exactly 7 wavelet scales 2 j Δt (j=0,1,…,6); the shortest scale is 2 0 Δt = Δt, the longest scale is 2 6 Δt = 64Δt = t MAX . The signal can be decomposed in terms of these dyadic time scales. We will define later (in sub-section 3.4.3) the concepts of wavelet "detail", "approximation" and "residual", at any given dyadic scale.

Continuous wavelets (a brief introduction)

A signal X(t) is represented formally via wavelet coefficients , X s C τ in continuous time as follows:

, , 2 0 1 1 ( ) ( ) X X t C t d d K σ τ σ τ σ τ ψ σ τ ψ τ σ σ =+∞ =+∞ = = -∞   =         (13a) , 1 ( ) t t σ τ τ ψ ψ σ σ -   =     ; t∈ℝ; τ∈ℝ; σ∈ℝ-{0 (13b) , , , ( ) ( )
( ), ( ) X C Xt td t Xt t σ τ σ τ σ τ ψ ψ +∞ * * -∞ = ≡  (13c)
where "t" is time, "τ" is a time-shift parameter, and "σ" is a scale parameter (a time scale parameter). It is assumed, for simplicity, that the signal is known for t∈[-∞,+∞]. Eq. 13a gives the decomposition of signal X(t). Eq. 13b expresses the (dilated) wavelet function at a given scale σ. Eq. 13c gives the wavelet coefficients of the signal: they are the scalar products of X(t) with the (conjugate) dilated and translated wavelet functions ψ σ,τ *(t).

A discrete time analogue of these equations can be constructed as follows: (a) replace "t" with t(i) = i×Δt (i = …,-2,-1,0,+1,+2,…) where "i" is the time index; (b) replace "τ" with τ(k)= k×Δt (k =…,-2,-1,0,+1,+2,…) where "k" is the time-shift index; and (c) replace "σ" with σ(s) = s×Δt (s = 0,1,2,3…) where "s" is the scale index (discrete but not dyadic --see further below).

Similarly, a finite time window analogue can also be constructed. Thus, one can deal with the practical case of discrete time step Δt and finite window [-T MAX /2, +T MAX /2] or [0, T MAX ].

Different wavelet bases can be obtained depending on the choice of the "mother wavelet" ψ(t). However, in the framework of "continuous" wavelets, it is not possible to find a set of wavelet functions ψ s,k (t) that form an orthogonal basis. For example, the Morlet wavelet basis is not orthogonal, implying that there are redundancies in the Morlet wavelet decomposition.

Multi-resolution wavelets and dyadic time scales

In contrast, with dyadic multi-resolution wavelets, an orthogonal wavelet basis with compact support can be constructed [START_REF] Daubechies | Orthonormal bases of compactly supported wavelets[END_REF][START_REF] Daubechies | Ten lectures on wavelets[END_REF], leading to multi-resolution wavelet analysis (below) 7 . This is made possible by the use of dyadic scales of the form σ(j) = 2 j ×Δt (powers of 2). Thus, in what follows, let "i" be the time index; "j" the dyadic scale exponent (power); and "k" the time-shift index. The decomposition of X(t) on the multiresolution wavelet basis takes the form:

, , 0 ( ) ( ) j k X j k j k j k X i C i ψ =+∞ =+∞ = =-∞ =  
(where "i" can be replaced by "t i ") (14a) 2 , ( ) 2

(2 )

j j j k i i k t t ψ ψ τ = × - (14b) , , , ( ) ( )
( ), ( ) i X j k i j k i i j k i i C Xt t Xt t ψ ψ =+∞ * * =-∞ = ≡  (14c)
where: j ≥ 0 : exponent of dyadic scale dilatation "2 j ו" (or compression if j ≤ 0); k: discrete time shift ↔ τ(k)= k×Δt; i: discrete time ↔ t(i)= i×Δt; i∈ℤ; k∈ℤ; j∈ℕ (the latter is the discrete dyadic scale)

The ψ j,k (t) basis functions form an orthonormal basis, with their images obtained by time translations and scale dilations/compressions of the mother wavelet ψ(t) 8 . Thus, we have:

7 Two functions f(t) and g(t) are "orthogonal" if their scalar product <f(t),g(t)> is null, where we define <f(t),g(t)> as the integral of f(t)×g(t) over the time domain. A function f(t) has compact support if it is null outside a finite interval [a,b]. 8 In the literature, the decomposition of X(t) is formulated in such a way that only dilatations (or "dilations") are retained, starting from the smallest scale Δt at j = 0, and increasing "j"; but it is possible to formulate this differently in terms of both dilatations and compressions.

m n m n m m n n t t d t ψ ψ δ δ +∞ -∞ =  ( 15
)
where δ is the classical Kronecker symbol defined by δ ij =1 if i=j and δ ij =0 if i≠j. Eq. 14a shows that the discrete time signal X(i) or X(t(i)), of finite energy, can be broken up into a linear combination of translations and dyadic dilatations of the basic functions with adequate coefficients C(j,k) given by Eq. 14c. Note that Eq. 14a can be viewed as the synthesis of a signal based on wavelet coefficients C(j,k). It can also be interpreted as a succession of approximations of the signal X(i) (in the least-squares sense) by a sequence {Xn(i): n=1,…..N} defined by: 1 , , 0 ( ) ( )

j n k X n j k j k j k X i C i ψ = -=+∞ = =-∞ =   ( ) ( ) ( ) as n n X i X i → → ∞. ( 16
)
The above equations constitute the conceptual basis of multi-resolution wavelet analysis.

For more details on the theory, see [START_REF] Mallat | A Theory For Multiresolution Signal Decomposition: The Wavelet Representation[END_REF]. Below, we focus on some useful consequences of the orthonormal wavelet decomposition above, introducing the wavelet approximation + residual, and interpreting the wavelet "approximation" as a filter of given dyadic scale.

Multiresolution wavelet decomposition: approximation A(t), residual R(t), and detail D(t) or component C(t)

In this subsection, we define the important concept of wavelet "detail" or "component", i.e., and we show how to extract from the signal X(t) a wavelet "component" with any given dyadic scale. It is important to note that only dyadic scales can be extracted; however, this is the price to pay for the other "good" properties (compact support and orthogonality of the decomposition). We also present the decomposition of signal X(t) into an "approximation" A(t) plus a "residual" R(t), both defined at any given cut-off scale (again, the cut-off scale can only be dyadic). Remarkably, the wavelet "approximation" can be viewed as a "filtered" version of the original signal: it contains only longer time scale components (longer than the given dyadic cut-off scale).

Orthogonal decomposition of X(t) -wavelet components, approximation, residual

Let us assume for simplicity that the (finite) discrete signal X(t i ) has a dyadic length (N=2 M ), the number of time steps is a power of 2 (dyadic). The multi-resolution wavelet decomposition above can be expressed as a succession of nested "approximations" A m , corresponding to increasing scales "m" (thus, "m" can designate a chosen dyadic time scale cut-off).

The difference between the signal X(t i ) and its approximation of order "m" is the "residual".

The difference between the approximations of order (m) and (m+1) is called the "detail" of order (m+1).

Intuitively, the approximation corresponds to the smoothed image of the signal, while the "detail" is a particular wavelet component that highlights the irregularities of the signal at a certain scale.

Starting from Eq. 14a, we obtain at each selected dyadic scale "m" the "wavelet detail" D m (t i ):

, , ( ) ( )

k X X m i m k m k i k D t C t ψ =+∞ =-∞ =  ⇔ , , ( ) , .
( )

m X m k m k k D i X i ψ ψ +∞ = -∞ =  (17a)
where "t(i)" can also be replaced by "i". The approximation A m of the signal X(t i ) at dyadic scale resolution "m" (cut-off scale "m"), and the residual R m at the same scale "m", are then given by: 1 ( ) ( )

X X m i j i j m R t D t ≤ - =  ; ( ) ( ) ( ) X X m i i m i A t X t R t = -  ( ) ( ) X X m i j i j m A t D t ≥ =  (17b)
Thus, at any dyadic scale "m", the signal is the sum of its approximation and residual:

( ) ( ) ( ) X X i m i m i X t A t R t = + (17c)
Using the latter equations to calculate the approximation of order (m-1) leads to:

1 1 ( ) ( ) ( ) X X X m i m i m i A t A t D t - - = + ⇔ 1 ( ) ( ) ( ) X X X m i m i m i A t A t D t + = - (17d)
As "m" increases, the A m 's become coarser and coarser approximations. The above equation shows that the "coarser" approximation A m+1 is obtained from A m by subtracting the detail D m . Note that the zero order residual R 0 is null by construction. The zero order approximation A 0 is identical with the original signal. Thus, we have, with our index conventions: 0 ( ) 0

X i R t = ; 1 0 ( ) ( ) X X i i R t D t = ; ...; ( ) ( ) ( ) X X M i i M i R t X t D t = - ; 1 ( ) ( ) X M i i R t X t + = (17e) ( ) 0 ( ) X i i A t X t = ;
( )

1 0 ( ) ( ) X X i i i A t X t D t = - ;...; ( ) ( ) X X M i M i A t D t = ( 17f
)
where M is the exponent of the largest dyadic scale. For example, M = 6 if the signal length N is in the range 64 ≤ N ≤ 127. The dyadic wavelet decomposition is performed without loss of information if N is a power of "2", such as N = 64; but if N = 127 the loss is maximal.

Theoretical interpretation of multiresolution wavelets (nested spaces and filters)

For completeness, note that a scale function ϕ(t) is also introduced. It can be shown that the "approximation" can be directly decomposed orthogonally over the set of dilated/translated scale functions ϕ j,k (t) (not shown here). The scale function ϕ(t) is also called the father wavelet, to be distinguished from the mother wavelet ψ(t) used above in the decomposition of Eq. 17a.

Theoretically, approximations A m and details D m are each decomposed into bases of ϕ's and ψ's, respectively, which form nested spaces called V m and W m , respectively. Considering a given scale "m", and recalling that m ≥ 1 ("m" increases for increasing dyadic scale 2 m ×Δt), it can be seen that the "detail" D m represents the wavelet component of scale "(m)", the residual R m incorporates the sum of all scales strictly smaller than "m", and the approximation A m incorporates all scales greater or equal to "m" (note: the approximation "m" includes the scale "m" itself).

The nested spaces V and W are related by V m+1 = V m -W m , as indicated by Eq. 17d. In fact, Eq. 17d describes a sequence of multi-resolution filters (the A m 's) and their residuals (the R m 's).

Remarks on multi-resolution approximation as a filter

The process of selecting a given wavelet scale "m" is similar to passing a "band pass" filter on the signal. The wavelet "approximation" A m can be viewed as a "low pass filter" of given dyadic scale. This filter is somewhat similar to a moving average filtering of the original signal. It was shown empirically in some cases that the two filters are quite close. See for instance [START_REF] Wang | Signal Processing of Water Level Fluctuations in a Sloping Sandy Beach Modeled in a Laboratory Wave Canal[END_REF] concerning water level signals in coastal & beach hydrodynamics. The advantage of the wavelet approximation as a "filter" is that it has the property of being orthogonal to its residual.

Final remarks on wavelets notations and terminology (pitfalls)

The presentation above is certainly too brief to be complete. Therefore, let us emphasize a few more points in order to avoid any possible confusion with multiresolution wavelets.

i. The wavelet "detail" D m (t i ) defined above can also be seen as a wavelet "component", to be denoted C m (t i ): this notation designates the wavelet component (C) at the m thdyadic scale (m). The sum of all components (over all scales "m") yields the signal itself. Note: from now on, we will prefer to use this notation, component C m (t i ), rather than D m (t i ). ii. For each given scale "m", the wavelet "component" C m (t i ) is a discrete time signal. This signal (component C(t)) should not be confused with the "wavelet coefficients C jk ", defined earlier, where "k" designates a time shift parameter. The "content" of signal X(t i ), at dyadic scale "j", is the sum over all time shifts "k" of the product Cjk×ψjk(ti)…The coefficient Cjk alone is not a time signal. iii. Multi-resolution wavelets can be presented differently in other texts, with other index conventions and notations. For example, the dyadic time scale labeled "m" can either increase or decrease with index "m" (here we have chosen that the scale increases with "m"). Also, some authors choose (1/scale) instead of (scale) in their visualisations of the wavelet decomposition in scale-time space.

Example application: Multiresolution wavelet analysis of a synthetic signal

Given the limitations of Fourier analysis, multiresolution analysis is particularly useful as a tool for analysis, synthesis, and/or compression, of complex non stationary signals. For this reason, we present briefly a synthetic example involving a signal that has a simple harmonic component added to a strongly non stationary mean or trend, as follows: ( )

( ) ( ) ( ) ( ) 1
0,1, 2,..., 31 i = ; 1 t Δ = (18)
This signal has an exactly dyadic length N = 32 = 2 5 (the dyadic exponent is M = 5).

Figure 8 (a,b,c,d) shows the signal X(t) and three different decompositions using the wavelet bases of Haar (Daubechies 2), Daubechies 4, and Daubechies 20. In each case, the six wavelet components C i (t) are shown (i = 0, 1, ..., 5). It can be seen that component C 2 (t) tends to capture the sinusoïd of period T=8. The corresponding dyadic scale is S j = 2 j Δt = 4 with j=2. Therefore, we see that the wavelet scale corresponds in this case to the half-period of the signal. On the other hand, the non stationary mean level is captured by component C 4 (t), with scale S j = 2 j Δt = 16, which corresponds to the half-duration of the non stationary signal. Note also that the components become smoother as the degree of the wavelet basis increases. The Haar wavelet is well suited to capture abrupt changes of mean level, but not for capturing smoother type of non stationarity.

(t) = 0, R 1 (t) = C 1 (t), .
.., and R 5 (t) = X(t)-m X = X(t)-A 5 (t) (note that R6(t) = X(t), not shown). (b): Plot of successive approximations A 0 (t) = X(t), A 1 (t) = X(t)-C 1 (t), ..., and A 5 (t) = X(t)..

Methods and techniques of wavelet analyses and cross-analyses

The above multi-resolution wavelet theory can be used to study different aspects of the signals, depending on the objectives:

• One can study the "detail" or (perhaps more interestingly) the "residual" of a signal X(t) at any selected dyadic scale "m". The analysis can be performed for several signals X(t), Y(t). One can compare for instance the details of air and pore pressure at scale 24h.

•

One can select dyadic components from one or several signals, and these can be studied and compared as a function of dyadic scale (2 m *Δt) across the whole spectrum of available scales, from finest scale Δt to largest dyadic scale (which is never smaller than half the record length).

Wavelets can also be used in different ways for cross-analyses of two signals like atmospheric and pore pressure. For instance:

• It is possible to develop a somewhat sophisticated approach for orthogonal crossanalyses of signals based on the theory of wavelet spectra and cross-spectra (see for instance Labat 2005 and references therein), but this theory will not be used here.

•

Simpler analyses can be developed by studying the cross-correlation between atmospheric pressure (Patm) and pore pressure (P). The wavelet component of interest is treated as a signal, and cross-component correlation analysis can be conducted at the chosen time scale: see Labat et al. (1999 a,b); [START_REF] Wang | Signal Processing of Water Level Fluctuations in a Sloping Sandy Beach Modeled in a Laboratory Wave Canal[END_REF]; Möri et al. (2011); for applications of cross-wavelet analyses in karst hydrology, beach hydrodynamics, and clay rock hydrogeology, respectively.

•

Yet another practical approach is to compare directly the amplitudes of the wavelet components of the two signals X(t) and Y(t) at the same dyadic scale "m".

The latter alternative may be viewed as a "simplified cross-wavelet approach". It is a way to combine two "simple wavelet analyses" to finally obtain a simplified but practical crossanalysis of the two signals at the selected scale. In this work, we will use this type of approach to identify semi-diurnal earth tide effects in pore pressure data, leading to an estimate of the elastic "Ss" as the ratio of tidal strain / pore pressure amplitudes in terms of semidiurnal wavelet components.9

In Sections 4 & 5, this simplified cross-wavelet approach will be applied to obtain a constant global "Ss" and also to quantify the effect of excavation on a more local estimate of Ss(t). To achieve this, we will need to characterize the time varying amplitude of a semi-diurnal wavelet component (time scale of moon induced earth tides). In his way, we will be able to "see" the effect of the moving excavation front on the semi-diurnal wavelet "signal".

For this task, we need a new tool that can evaluate statistically the "time-modulated amplitude" of a signal. This tool is available already (although we apply it in an unusual way, in combination with wavelets): it is the theory of "statistical envelopes", presented in the next sub-section below.

Statistical envelope analysis of a signal ("modulated amplitude")

The Cramer-Leadbetter theory 10 and the related Hilbert Transform 11 , are the theoretical tools that enable us to define, mathematically and statistically, the "modulated amplitude" (or the "envelope") of a fluctuating and evolutionary pressure process. The theory of statistical envelopes is well established for zero mean stationary random processes (Vanmarcke 1983;[START_REF] Veneziano | Envelopes of Vector Random processes and their Crossing Rates[END_REF]. This theory, together with the associated Hilbert Transform, is particularly useful for characterizing relatively narrow-band signals, e.g. ocean & coastal waves [START_REF] Huang | A new view of nonlinear water waves: The Hilbert spectrum[END_REF].

But the present study (see Section 5) suggests that the "statistical envelope" can also be useful in a variety of other cases in geophysics, provided it is combined with carefully selected "components" of the signals (e.g., here, the semi-diurnal wavelet component of pore pressure). In fact, the mathematical theory of statistical envelopes is established only for stationary processes, but as will be seen in Section 5, using wavelets, we will get rid of the strong non stationary mean trend due to excavation. This is naturally accomplished by using multi-resolution wavelets at specific dyadic scales. The dominant scale of the mean trend is ∼11days, while the wavelet scale selected for analysis is semi-diurnal (scale ≈ 12hours). The effects of the 11 day trend due to excavation are still "visible" as a time varying "modulated" amplitude in the semi-diurnal wavelet component C 12HOURS (t), which we capture using envelope theory. Thus, our use of the theory of envelopes and Hilbert Transform in the present context, with semi-diurnal pressure fluctuations superimposed on a strong nonlinear trend due to excavation, seems rather new in the literature.

Cramer-Leadbetter theory in a nutshell.

For a zero mean (detrended) stationary random process X(t), one defines first a dephased "sister" process S(t). Based on X(t) and S(t), one obtains a symmetric double envelope ±R(t) such that:

R(t)=(X 2 (t)+S 2 (t)) 1/2 .
The envelope is usually much smoother than the original process X(t).

In the special case where X(t)=sin(ωt), we obtain S(t)=cos(ωt), and then R(t) = sin 2 (ωt)+cos 2 (ωt) = 1, as one may expect. Indeed, it should be obvious that the envelope of the sinusoïdal signal X(t)=sin(ωt) is just a pair of horizontal lines of ordinates -1 and +1.

Hilbert transform in a nutshell.

The Hilbert transform is an integral transform of the form:

( )

1 ( ) ˆ( ) ( ( )) ( ) ( ) ( ) X t X t H X t h X t X t h t d d t τ τ τ π τ +∞ +∞ -∞ -∞ = = ⊗ = - = -   .
For a real-valued signal X(t), H(X(t) is essentially the convolution of the signal X(t) by the kernel, or transfer function, h(t) = 1/(πt). Note:

X(t) = sin(t)  H(X(t)) = -cos(t).
The Hilbert transform can be implemented in a complex framework as follows (as in the Matlab software): let X(t) the real valued process; let Y(t) = H(X(t) its real Hilbert transform. Then the complex Hilbert transform of signal X(t) is defined as Z(t) = X(t)+i.Y(t) where i is the pure imaginary number (i = √-1). The above defined envelope of X(t) is then obtained as:

( ) ( ) ( ) 2 2 2 2 ( ) ( ) ( ) ( ) R t X t H X t X t Y t Z t = ± + = ± + = ± .
The Hilbert transform is implemented numerically for discrete time signals in MATLAB.

Box 1. Envelope theory in a nutshell

Rather than present here the full theory of statistical envelopes, we refer the reader to the abovecited references for the theory, and to the summary presented here in Box N°1. The reader who is not directly interested in the theory can move on to Section 5 for application of statistical envelope + wavelet analysis, to the non stationary evolution of semi-diurnal fluctuations of (relative) pore pressure during the excavation of Mont Terri gallery Ga98 ("mine by test").

As mentioned earlier, the final objective and result of this analysis is to estimate the evolution of specific storativity "Ss(t)" in the EdZ of the clay rock. This is done by using the tidal dilatation model of Bredehoeft relating tidal strain / relative pore pressure. The originality of this "envelope" approach should be emphasized. The amplitude of the semidiurnal wavelet component is estimated as a modulated envelope, which varies more smoothly than the signal itself. The resulting storativity coefficient is therefore a timevarying coefficient, Ss(t), which reflects the mechanical disturbance due to the evolving excavation front, over short time scales (days, weeks).

Physically-based hydro-mechanical models to characterize the clay rock from pressure signals (earth tides and barometric effect)

The statistical methods of signal analysis are complemented or combined with physically based, but "simplified" hydro-mechanical models. These are used for hydrogeologic interpretation of the pressure signals (pore pressure and atmospheric pressure), and they are also related to inferred earth tide strain fluctuations, with the aim of characterizing quantitatively some of the hydrogeologic properties of the porous clay rock (Ss, Φ).

Specifically, the simplified hydro-mechanical models used here relay on assuming both water and the prous matrix (bulk) to be elastically compressibility, and on accepting Terzaghi's 1936 effective stress concept. The corresponding models can be traced back to [START_REF] Jacob | On the flow of water in an artesian aquifer[END_REF] for the barometric efficiency model; and to [START_REF] Bredehoeft | Response of Well-aquifer Systems to Earth Tides[END_REF] for the earth tides model (also [START_REF] Jacob | On the flow of water in an artesian aquifer[END_REF]Cooper 1966). These physically-based models both assume that the total stress can be decomposed additively into effective stress between grains, and pressure (fluid stress), according to the Terzaghi concept. This implies that the grains of the claystone are "incompressible". The more general Biot poroelasticity theory can account for the case where individual grains are compressible (see for instance Hsieh et al. 1988 concerning specific storativity).

In this work we use Terzaghi based hydro-mechanical models. They are outlined below in Box N°2, where the final results are expressed in terms of:

• The specific (elastic) storativity Ss [m -1 ]; •
The effective (elastic or dynamic) porosity Φ [m 3 /m 3 ].

Long time scales: statistical structure of pressure signals and characterization of time averaged properties (specific storativity Ss; elastic porosity φ)

This section is devoted to the analysis of long records (months to years) of observed air pressure and relative pore pressure signals, for the purpose of characterizing global and average hydro-geologic properties of the clay stone (global in space, and also, average in time)…Statistical analyses based on correlations, spectra, and wavelets, will enable us to detect and characterize in particular the diurnal relation between relative pore pressure and air pressure (leading to a "barometric efficiency"), and the half-diurnal effect of earth tides (leading to a characterization of compressibility and/or specific storativity).

The purpose of this "box" is to summarize the simplified, but physically-based, hydro-mechanical models. These models are used in combination with statistical analyses in order to interpret pressure signals in terms of rock properties (elastic storativity Ss, effective elastic porosity Φ).

Tidal strain -astronomic and global earth (the semi-diurnal component)

Astronomic movements of Earth, Moon, Sun, induce fluctuations of the gravimetric field g(t), which has an impact on deformations of the earth crust (displacements, strains). Focusing on the M2 harmonic (Moon semi-diurnal), the volumetric strain fluctuations take the form:

( ) ( ) 2 .sin M t t ε ε ω φ = Δ + where ε Δ is the tidal strain amplitude [m 3 /m 3 ].
Concerning M2, it is the principal lunar earth tide, and its period is more precisely 12h 25 mn. [START_REF] Takeuchi | On the earth tide of the compressible earth of variable density and elasticity[END_REF] gives:

( )

2 0.5 W ag ε Δ = ×
, where W2 is the tidal potential of M2, "a" is the radius of the earth, and "g" is gravity. One may take W2/g ≈ 10 for M2. The constant 0.5 was computed by Takeuchi from earth crust mechanics data (Lamé stiffness moduli, Love number). At the scale of a geologic unit (as opposed to the global scale of the earth), [START_REF] Melchior | The tides of the planet earth[END_REF]

estimates 8 3 3 2 10 [ / ] m m ε - Δ ≈ ×
(value adopted in this work in the absence of direct measurements).

Tidal strain and elastic specific storativity Ss of a geologic unit -the Bredehoeft model

Specific storativity Ss : Ss is in units of [m -1 ] ↔ Ss/(ρg) is in units of [(m3/m3)/Pa].

( )

1 [ ] S M P M P S m V V H δ δ -=
where VMP is the volume of deformable bulk Porous Matrix, and δH is total head variation, with head H defined by:

ATM w p p H z g ρ - = + . Equivalently: ( ) ( ) ( ) 1 [ ] S M PM P R E L S g Pa V V p t ρ δ δ -= , with ( ) ( ) ( ) REL ATM p t p t p t = - (relative pressure).
Using Terzaghi's 1936 effective stress concept, assuming constant total stress, and assuming elastic compressibility of water (β=1/Kw) and of the bulk porous matrix (α=1/Kb…) leads to:

( )

1 S w w B W S g g K K ρ α β ρ   Φ = +Φ = +    
(where Kw 2 GigaPa for water).

The [START_REF] Bredehoeft | Response of Well-aquifer Systems to Earth Tides[END_REF] .

Barometric efficiency B -classical model (compressibility & effective stress)

Without going into details (see [START_REF] Freeze | Groundwater[END_REF], we present the result obtained with the effective stress hypothesis [START_REF] Terzaghi | The Shearing Resistance of Saturated Soils and the Angle between the Planes of Shear[END_REF], in the case where the total stress is held constant, and assuming here also that our hydro-geologic unit behaves like a (partially) confined aquifer:

( ) ( )

REL ATM p t B p t δ δ = -×
, with the barometric efficiency B given by: ( )

B β α β = Φ + Φ .
Using the previous expression for Ss in terms of (α,β) compressibilities

( ) ( ) s B S g β ρ  = Φ .
In this work, barometric efficiency is re-interpreted statistically as the spectral gain at frequency 1/24h. The above equation shows that barometric efficiency B, specific elastic storativity Ss, and effective elastic porosity Φ are interrelated (the compressibility β=1/Kw of water is known).

Box 2. Summary of physically based hydro-mechanical models

Long records with many time steps (tens of thousands or much more) usually have non stationary structures at various scales. In this section, the undesired non stationary structures are partially eliminated by the selection of specific peak frequencies (in Fourier spectra), and/or by the selection of pertinent dyadic "scales" via multi-resolution wavelet analysis. We seek here only global characteristics of the signals (global amplitudes, etc.).

As explained earlier in the methodology Section 3, in order to achieve these goals (hydrogeologic characterization), we use simplified hydro-mechanical "models" to account for the porosity "Φ", the compressibility of water, and the compressibility of the bulk porous matrix. This leads also to the related concept of elastic specific storativity "Ss", based on Terzaghi's decomposition of total stress into effective stress+pressure. The "Ss" coefficient represents an equivalent elastic coefficient for the water filled poro-elastic rock. See the "outlook" comments in Section 6 for a discussion on the limitations and possible extensions of these hydro-mechanical models.

A first example: statistical analyses and interpretations of "clean" records of raw (un-processed) pressure signals over∼1 month

The above-mentioned statistical techniques are first illustrated here using only relatively short sections of the pore pressure signal PP1(t) at Mont Terri, with constant time step and no data gaps (no reconstitution required): see Figures 10, 11, 12, 13. The unprocessed pore pressure signal is analysed in terms of:

1. its frequency spectrum (shown in Figure 11), and 2. its wavelet component at nearly semi-diurnal time scale (shown in Figure 13).

Specific storativity from earth tides (using semi-diurnal wavelets)

We estimate in this subsection the month-scale specific storativity "Ss" from semi-diurnal wavelet component of relative pore pressure (strain/pressure --earth tide effect).

The specific storage coefficient quantifies the volume of water stored per unit volume of porous domain and per unit variation of hydraulic head. It can be estimated from earth tide amplitude using the relation of [START_REF] Bredehoeft | Response of Well-aquifer Systems to Earth Tides[END_REF], based on a simplified model that takes into account bulk deformation effects due to earth tides in the poro-elastic water filled porous medium:

s S h ε Δ = Δ (19)
where S S is specific storativity (m -1 ), |Δε| is the amplitude of volumetric strain fluctuations related to the M2 semi-diurnal earth tide, and |Δh| is the amplitude of relative pressure head fluctuations (estimated from the semi-diurnal wavelet component of relative pressure). The strain amplitude |Δε| was estimated by [START_REF] Takeuchi | On the earth tide of the compressible earth of variable density and elasticity[END_REF] based on astronomic earth tides and earth elasticity, then later on (in the 1960's) by Melchior and others for hydro-geologic formations: see [START_REF] Melchior | The tides of the planet earth[END_REF]Melchior ( ) (originally in 1960)), Bredhoeft (1967), and others: |Δε| = 2×10 -8 m 3 /m 3 .

Note that the site specific signal of volumetric strain Δε(t) is not available, so its amplitude is only estimated indirectly (from data surveys in the literature). However, our method can be extended to the case where Δε(t) is available, by applying a semi-diurnal wavelet analysis on the Δε(t) signal.

The amplitude of relative pressure head Δh(t) is obtained from multi-resolution wavelet analysis of P(t)-P ATM (t). The wavelet component with dyadic time scale close to semi-diurnal is isolated, and its amplitude is evaluated using either of the following "norms":

• the RMS (Root-Mean-Square) standard deviation: "norm 2" or 2 ; or else • the absolute mean deviation ("norm 1" or 1 ). Note:

1 ( ) 1 ( ) p p p X i X N σ   = -       (1) 1 ( ) X i X N σ = -  and ( ) 2 (2) 1 ( ) X i X N σ = -  (20)
Figure 12 shows the evolution of the components according to time and scale. The goal is to isolate the component corresponding to earth tides (the semi-diurnal component): as indicated on this figure, we have chosen the component of dyadic scale T = 8h (close to 12h). Note: in other calculations, we have attempted to improve the accuracy of semi-diurnal scale selection by computing the mean of the dyadic 8h and 16h wavelet components: Figure 13 shows the evolution of the selected wavelet component (in this case, C 8H (t)). We have determined numerically the amplitude of the fluctuations of this wavelet component, e.g., its absolute mean deviation. The result can be expressed in terms of relative pressure head: for instance it was found that the absolute mean amplitude ("norm1") of semi-diurnal pressure head is Δh = 0.058 m. We can then calculate specific storativity Ss by applying the previous relation (Eq. 19). The results are the following, using either "norm1" or "norm2" for Δh:

(C 8H (t)+C 16H (t))/2.
Δh(norm 1)=0.058 m; S s (norm 1)= 3.4×10 -7 m -1 Δh(norm 2)=0.064 m ; S s (norm 2)=3.1×10 -7 m -1

Effective elastic porosity Φ from barometric efficiency (using spectral gain)

We now estimate month-scale effective porosity from the diurnal frequency gain between relative pore pressure and air pressure (barometric efficiency).

The diurnal (and other) components would be due in part to barometric pressure effects. The effect of barometric pressure fluctuations on observation wells in closed, confined aquifers, is well known. The piezometric head in the aquifer is related linearly to atmospheric pressure (assuming linear poro-elastic behavior) via a coefficient named "barometric efficiency" (B or B E ).

Although the geometry, boundary conditions, and hydraulic configuration of the Mont Terri claystone appear quite different from those of a confined aquifer, we nevertheless propose to retain the same linear relationship. This makes it possible to calculate effective porosity, according to the relations described in Section 3.6 and Box N°2 (the following is Jacob's 1940 relation):

w S E B S g ρ × Φ = × × (21)
where Φ or Φ EFF is the effective elastic porosity (m 3 /m 3 ), B the barometric efficiency which expresses the elastic response of the system, ρ the density of water, g the acceleration of gravity, E W (or K W ) is the stiffness modulus of water (E W = K W ≈ 2.05 GigaPascals).

In order to calculate barometric efficiency B, a statistical cross-spectral analysis of the relation between atmospheric pressure and relative pressure was developed, with 30 minute time steps, to determine the reality of the effect of barometric pressure. We carried out the cross-spectral analyses between barometric pressure and relative pressure in borehole BPP-1 (from the PP niche, which contains two rooms, PP1 and PP2).

The estimated spectral gain function (in frequency space) was finally used to determine barometric efficiency (B) by looking at the value of the spectral gain at diurnal frequency, as shown in Figure 14. From this spectral gain, barometric efficiency is determined, and porosity is calculated with Eq. 21 just above, using the previous estimate of specific storativity (Ss) from semi-diurnal wavelet analysis. We obtain in this fashion the following results:

φ EFFECTIVE (norm 1) ≈ 6%; φ EFFECTIVE (norm 2) ≈ 5.4%
Recall that this procedure to obtain "Φ" (by spectral analysis of barometric gain at diurnal frequency) is combined with the previous procedure to estimate "Ss" (by wavelet analysis of semi-diurnal earth tide effects). In fact, the two procedures do not need to be implemented sequentially. They can be viewed as a system of two equations (Eq. 19 and Eq. 21) with two unknowns (Ss, Φ). The two unknowns could have been determined both at once (alternative approach).

Long term analysis of pressure signals over more than 1 year (after pre-processing)

The same techniques illustrated in the previous Section 4.1 are now being applied to a much longer -but pre-processed -joint record of absolute pore pressure and atmospheric pressure (p(t), p ATM (t)). The time span of the processed pair of signals is about 15 months (from 29 January 2004 to 12 April 2005).

Note that, initially, only shorter pieces of "clean" signals were available during that time window. Note also that the pre-processing procedures (reconstructing, homogenizing, etc.) must be implemented jointly on p(t) and on p ATM (t) over the same time window, since the relative pressure signal must be obtained from p REL (t) = p(t)-p ATM (t), and cross-analysis of (p REL (t), p ATM (t)) must be implemented to identify the barometric effect. A flow-chart depicting pre-processing tasks and methods was shown earlier in Sub-Section 2.2.1 (Figure 2).

The results of statistical analyses of the 15 month long pre-processed pressure signals (p REL (t), p ATM (t)) are shown below in Figures 15,16,17, 18 for the PP1 pore pressure sensor.

Figure 15 shows p REL (t) over the entire 15 month time span. Similarly to the previous Section 4.1, we analysed the cross-spectral gain of (p REL (t), p ATM (t)) at frequency 1/24h (not shown here), and the wavelet component at nearly semi-diurnal time scale (Figure 17). In addition, Figure 16 shows the spectrum of Δp REL (t), the first order difference of p REL (t). Differencing enhances the peak frequencies, and allows one to detect a dominant 1/24h frequency, and a somewhat less dominant 1/12h frequency (subdominant but still visible). We see that these two characteristic frequencies can be detected in spite of the obvious non stationarity (large scale non linear trend) of the 15 month pore pressure signal (Figure 15). Φ (norm 1) ≈ 1.0% |Δh| (norm 2) ≈ 0.08m ; Ss (norm 2) ≈ 2.6×10 -7 m -1 Φ (norm 2) ≈ 0.9%

where "norm 1" stands for the absolute mean, and "norm 2" stands for the Root Mean Square (or Standard Deviation). We focus here on a total duration of ∼5 months, corresponding to the syn-excavation period of gallery Ga98 (17th of November 1997 to 23rd of April 1998). The pre-and post-excavation periods were also analysed for comparison [START_REF] Fatmi | Méthodologie d'analyse des signaux et caractérisation hydrogéologique: Application aux chroniques de données obtenues aux laboratoires souterrains du Mont Terri[END_REF], but it is sufficient to treat here the 5 month "syn-excavation" period. Furthermore, we will also focus, within this 5 month period, on a shorter phase of about 2 weeks of "fast" pressure changes (see below). 19, it is clear that the pore pressure signal in sensor PP1 exhibits a strong response to the passage of the excavation front. This indicates the existence of a moving EdZ ("Excavation disturbed Zone") as the excavation front passes near the pressure sensors. Pressure rises gradually at first, then very rapidly ("fast" 11 day period). After this rise, pressure seems to reach a stable value but in fact, it starts a long and slow relaxation process (confirmed from analyses of longer records). The sharp rise of pore pressure can be explained as a hydro-mechanical disturbance of the geologic porous formation due to excavation works. More precisely, Figure 19 shows the pore pressure data from the PP1 sensor, in comparison with the evolving distance D(t) of the excavation front. The pre-and post-excavation phases are not shown. The 5 month "syn-excavation" phase was determined to run from the 17th of November 1997 to the 23rd of April 1998, and it is shown entirely in this figure (see caption for other details and comments). Clearly, within the 5 month syn-excavation phase, there is a remarkable sub-phase of "fast" pressure rise, lasting about two weeks only (∼11 days). Similar effects were observed for the other sensor PP2 in the same borehole (not shown here). A study of the signal on longer time scales (not shown in this text) confirms in fact that the pore pressure disturbance was indeed reversible (hence the name "Excavation disturbed Zone", as explained earlier). Indeed, a slow relaxation of pore pressure was observed during a long period, somewhat less than two years, subsequent to the end of excavations works (there are almost two years counting from the start of excavation till the end of pressure relaxation: 30/11/1997-18/02/2000). The pore pressure at the end of this period reached approximately the value it had at the beginning of the period (pre-excavation) and it continued to decrease to reach finally the pre-pumping test value (see Fatmi 2009, Fig.28, p.93, for more details).

Short time scale effects of gallery excavation on

Upon inspection of Figure

These observations justify our initial assumption and remark on the distinction between "EDZ" (irreversibly damaged "Excavation Damaged Zone"), and "EdZ" (Excavation disturbed Zone): the latter is thought to be a reversibly disturbed zone due to quasi-elastic hydromechanical variations of pressure and stress. These reversible variations occurred in the zone of our pressure sensor, and they caused the observed variations of pressure during the passage of the excavation front (which came within 4 m of the pressure sensor at the closest point).

Now, looking at the short and fast rise period of 11 days in Figure 19, it may seem "impossible" at first sight to detect a pattern of semi-diurnal pressure fluctuations due to earth tides, and even less, to follow possible changes in these fluctuations during the phase of rapid pressure change. Yet, although earth tides are partially masked, they can still be detected indirectly using proper filtering techniques. For example, applying differential filters can enhance subdominant frequencies, moving average filters can be used to look at residuals, and multi-resolution wavelet analysis can be used to extract precisely the semidiurnal earth tide component (as will be seen).

Below, we will take advantage of the precise way in which multi-resolution wavelets can capture the semi-diurnal component of the original signal, and we will then analyse the new "12h" signal using envelope techniques, to finally obtain time-varying hydraulic parameters of the claystone during the passage of the excavation front nearby the pore pressure sensor.

Multi-resolution wavelet analysis of relative pressure signals during excavation: extraction of earth tide semi-diurnal component

First, we obtain the relative pore pressure p REL (t) from the absolute pore pressure signal p(t) [shown in Fig. 19 above as "PP1(t)"] by subtracting the atmospheric pressure signal p ATM (t) (not shown), whence: p REL (t) = p(t)-p ATM (t).

Obviously, one should not expect p ATM (t) to be influenced by excavation works (rather, p ATM (t) contains barometric fluctuations, with peak frequency 1/24h, part of these having an impact on pore pressure).

The procedure to identify the elastic specific storativity Ss is then similar to that illustrated in Section 4.1, except that here, after extracting the semi-diurnal wavelet component "C 12H (t)" from the relative pore pressure signal p REL (t), we perform an additional analysis of the modulated amplitude, or "envelope" of C 12H (t), in a statistical framework.

Note, from now on, the semi-diurnal wavelet component "C 12H (t)" is treated as a signal. The wavelet signal C 12H (t) and its statistical envelope are now presented and analysed in the next section below (Section 5.3), in order to derive, from that information, the time-varying specific storativity "Ss(t)" during the syn-excavation phase.

Statistical envelope of semi-diurnal wavelet of relative pressure during excavation (Hilbert transform & method of maxima)

In this section, we present the wavelet component signal "C 12H (t)" of relative pore pressure, and we analyse its modulated amplitude by computing its Cramer-Leadbetter statistical envelope using MATLAB's Hilbert Transform (Section 3.5: Statistical envelope analysis...).

In addition, in some cases, we compare this theoretical envelope with an empirical one computed by our own "method of maxima". Briefly, the latter method generates an empirical envelope that has irregular time steps. The result of the empirical method agrees generally with the theoretical envelope except for isolated discrepancies (this is a minor point, but the interested reader can consult Fatmi 2009 for more details).

Figure 20 shows (in blue) the semi-diurnal signal "C 12H (t)", i.e. the 12h wavelet component of relative pressure in sensor PP1, and it shows also (superimposed in red) the statistical Cramer-Leadbetter envelope of "C 12H (t)" (envelope computed with the Hilbert Transform).

Recall our gal is to analyse the effects of excavation of gallery Ga1998. Here, the entire "synexcavation" phase of 5 months is shown (approximately from mid-november 1997 to midapril 1998).

Figure 21 shows a temporal zoom of Figure 20 over a much shorter period. The time window is approximately one month, and it is centered around the "fast" sub-phase of synexcavation, which lasts only 11 days. Comparing Figure 21 to Figure 19, it is now clear (and remarkable) that the 11 day phase of rapidly increasing pore pressure (Figure 19) corresponds to a significant increase of the amplitude of the semi-diurnal component of pore pressure (Figure 21). Recall, from Fig. 19, that this 11 day phase corresponds to the approach of the excavation front nearest to the PP1 pressure sensor (distance D(t) ∼ 4 m at the closest point).

Calculation of specific storage Ss(t) as a function of time based on Bredehoeft's simplified hydro-mechanical model (strain / pressure)

By analogy with the procedure described in Section 4.1-4.1.1, we estimate again the elastic specific storativity Ss from the semi-diurnal wavelet component C 12H (t) of relative pore pressure (using now the average of the 2 components at dyadic scales 8h and 16h, respectively). And again, we invoke the Bredehoeft earth tide strain/pressure model (see Eq. 19 earlier).

Only, this time, instead of taking a constant amplitude for C 12H (t), we use instead the time-modulated amplitude of C 12H (t) represented by its statistical envelope "R(t)". On the other hand, since we still do not have direct in situ measurements of the volumetric strain Δε(t), we must replace Δε(t) by an indirect estimate of its amplitude based on data surveys: accordingly, we use again |Δε| = 2×10 -8 m 3 /m 3 (as in Sections 4.1 and 4.2). In this way, we obtain, instead of Eq. 19:

( ) ( )

S S t R t ε Δ = (19b)
Figure 22 shows the time evolution of log 10 Ss(t), where Ss(t) was estimated using the empirical envelope R(t) computed by the method of maxima, rather than the theoretical Cramer-Leadbetter envelope computed by the Hilbert Transform. The evolution of the excavation front distance D(t) is also shown on the same graph.

Figure 23 shows the time evolution of log 10 Ss(t), using the theoretical Cramer-Leadbetter envelope (computed by the Hilbert Transform) to estimate Ss(t). The evolution of the excavation front distance D(t) is also shown on the same graph.

Comparing Figure 22 and Figure 23, it is somewhat re-assuring that the two methods of estimating the wavelet envelope yield approximately the same results. Thus, we focus mostly on results obtained with the Hilbert Transform method for the remaining plots and analyses.

The time scale of Figure 23 is about 1.5 month, including the interesting 11 day window of fast changes. It can be seen that the specific storativity Ss(t) decreases by one to two orders of magnitude during just a few days, as the excavation front approaches within roughly 6 m to the PP1 pressure sensor (before reaching an even closer distance of 4 m).

Figure 24 shows a zoomed view of Ss(t) with a cartesian ordinate rather than semi-log. The zoomed time window is only about 2 weeks or so, to see more details in the 11 day period of fast change.

Finally, Figure 25 gives a broader picture, depicting the time evolution of log 10 Ss(t) over the entire "syn-excavation" phase of 5 months (including, again, the 11 day period of fast changes).

Syn-excavation Ga98

Fig. 24. Zoomed view of the evolution of Ss(t) on a cartesian ordinate (rather than logarithmic). The plot is centered on the 11 day period of fast pressure change, corresponding to the approach of the excavation front of gallery Ga98 to the pressure sensor. Finally, from the previously estimated evolution of Ss(t), and using the spectral barometric gain at 1/24h frequency (as in Section 4.1.2), it is possible to apply again the simplified hydro-mechanical model of Eq. 21, in order to estimate the time-varying elastic porosity Φ(t):

( ) ( )

w S E B t S t g ρ × Φ = × × (21b)
where Ss(t) is the time varying specific storativity obtained earlier from the strain/pressure envelope (Eq.11b, Section 5.4).

It should be observed that, in this analysis, the elastic porosity Φ(t) varies in time proportionally to Ss(t). It can be seen that a constant barometric efficiency "B" is assumed; the diurnal gain "g" between atmospheric and relative pore pressure is assumed constant during excavation. But, because pore pressure is influenced by excavation (even if air pressure is not), the assumption of constant barometric gain ("B" or "g") during excavation is debatable and should be verified.

In conclusion, a revised version of the procedure should be to analyse the diurnal barometric gain, like the semi-diurnal earth tide effects, in a wavelet framework rather than in a Fourier spectral framework. This would provide a tool to either capture or eliminate, depending on purpose, any non stationary / evolutionary phenomena.

Summary, conclusions, and outlook

Summary and conclusions

We have presented in this chapter a set of statistical methods for pre-processing and analysing multivariate hydro-geo-meteorological time series, and we have developed hydro-mechanical interpretations, in order to infer hydro-geological properties such as Ss and Φ. The aim of this study is to quantify some of the properties of the clay rock in terms of elastic response to barometric and earth tides fluctuations, both globally and locally in spacetime. We focused on atmospheric and pore water pressure signals measured in Opalinus clay stone, at the Mont Terri underground rock laboratory (URL), in the framework of Mont Terri's LP 14 experiment.

Pre-processing methods were developed in order to detect and correct some defects common to most time series data, like aberrations or outliers, variable time steps, and missing data. These pre-processing steps are necessary, in most cases, to obtain signals that are both "clean" and long enough for useful analyses. Durations of a couple of weeks to a few months are needed for short scale analyses related to excavation works. Ideally, several years are required for global characterizations, and even more for climate trend analyses.

Once the pressure signals were pre-processed (reconstructed), statistical analyses were performed with cross-correlation functions and cross-spectra, and with multi-resolution wavelet analyses. These techniques were then used to determine "globally" some hydromechanical parameters of the Opalinus clay: elastic specific storativity (Ss), and effective elastic porosity (Φ). The values compare favourably with data from the literature, obtained independently with completely different methods. In particular, the simplified hydromechanical model of Bredehoeft relating p REL (t) to earth tides was applied to the "12h" dyadic component of p REL (t), leading to a global estimated Ss value of 10 -6 m -1 (1E-6 m -1 ) in excellent agreement with that deduced from hydraulic tests at the Mont Terri site.

After this "global" characterization of the clay from long term pore pressure measurements, we turned our attention to the (significant) impact of excavation works on pore pressure signals. We devised an original method to infer the time evolution of clay rock properties (Ss(t)) due to the excavation disturbance. The pressure disturbance coïncides (with some lead time or negative delay) with the arrival near the pressure sensor of the excavation front D(t) (or of the "EdZ" around and ahead of the front).

Indeed, pressure signals recorded in the PP1 chamber of the BPP-1 borehole during excavation of the Ga98 gallery, have shown a net increase of absolute pressure (+5.5 bars) over the 11 day period of fast change (or as much as +7 bars over the total 5 month "synexcavation" phase). This is evidence of a direct hydro-mechanical disturbance of the rock during the passage of the excavation front.

A study of the signal on longer time scales (from pre-to post-excavation stages, not shown in this text for lack of space) confirmed that the disturbance due to excavation was indeed reversible. A slow relaxation of pore pressure was observed during a period of about two years (or somewhat less) subsequent to the end of excavation, reaching approximately the value it had before excavation. This observation appears to justify our initial remark on the distinction between "EDZ" (irreversibly damaged "Excavation Damaged Zone"), and "EdZ" (Excavation disturbed Zone). The latter "EdZ" is a reversibly disturbed zone due to quasi-elastic hydro-mechanical variations of pressure and stress.

In summary, the signal analysis and processing methods applied at various (long/short) time scales, exploit the influence of natural "forcings" like earth tides and barometric fluctuations on pore pressure signals. Combining multi-resolution wavelet and he Cramer/Hilbert envelope analysis reveals a time modulation (disturbance) of the "12h" dyadic wavelet component, suggesting a modification in time of the specific storativity Ss(t) during the excavation of gallery passing nearby the pressure sensor. During the passage of excavation, there is a net decrease of Ss(t) by almost two orders of magnitude, from 10 -6 (1E-6) m -1 , down to 10 -8 (1E-8) m -1 . This minimum occurs during a short period of just a few days, corresponding to front crossing (actually, in advance of front crossing by a couple of meters).

This result demonstrates how statistical analyses on pore pressures can be a valuable tool in the assessment of potential clay repositories for the purpose of high level radioactive waste disposal in deep, isolating geologic formations.

Outlook

A useful extension of this work would involve the development of a multivariate or "multi-cross analysis" of time series data {Xi(t), i = 1 , 2, …, K}, where the Xi's are either reconstructed or raw signals, involving possibly a number of pore pressures measured in different boreholes at different measurement sections, as well as other types of signals such as air pressure, air temperature, air humidity (relative or absolute), wall crack displacements, and moisture contents in the rock.

Accordingly, in multi-cross analysis, the results would be expressed in matrix form M XiXj (s) where "s" is the relevant parameter (time lag, frequency, etc.): thus M XiXj (s) may represent a K×K matrix of cross-correlation functions, or of Fourier cross-spectra, and/or, of multiresolution cross-wavelets (or their cross-correlations at given scales). These analyses will be useful for assessing the added value of all signals collected at the site, in terms of the "information" gained concerning the hydro-geologic properties of clay rock at different time scales (hourly, daily, yearly).

The purpose of such multi-cross analyses might be to integrate the complete set of all available signals in order to solve an inverse problem, namely, the identification of hydrogeologic properties from observations of state variables. Thus, in the work presented here, the state variables were the pore pressure at different sensors, atmospheric pressure, and (indirectly) tidal strain. The properties to be identified were (Ss, Φ) -but this set of unknown coefficients could be extended in future, in a generalized treatment of the "hydro-geologic identification problem".

For instance, techniques to estimate the hydraulic conductivity "K" from natural signals sampled at different positions might be developed as a complement to classical well pumping tests (e.g. pulse tests). Indeed "K" is an important property for a geologic radioactive waste repository, as "K" must be low to efficiently isolate radioactive waste from the environment. Furthermore, it is well known that "K" may vary greatly in geologic media, depending on the spatial scale of analysis (e.g. [START_REF] Ababou | Quantitative Stochastic Hydrogeology: the Heterogeneous Environment. Chap. 8 in Part III of "Overexploitation & Contamination of Shared Groundwater Resources: Management, (Bio)Technological, and Political Approaches to Avoid Conflicts[END_REF][START_REF] Sahimi | Flow phenomena in rocks: from continuum models to fractals, percolation, cellular automata, and simulated annealing[END_REF][START_REF] Gelhar | Stochastic Subsurface Hydrology[END_REF][START_REF] Ababou | Self-Similar Randomness and Spectral Conditioning: Analysis of Scale Effects in Subsurface Hydrology. Chap. XIV[END_REF]. With the "multi-cross" concept, the spatial scales of analysis can be associated with the distances between sensors.

The multi-cross analysis of these signals will also serve as a complement to physically based models, with some of the signals considered as inputs, and others considered as outputs of the models. It will be very useful to link the statistical approaches with physically-based models that can account for the hydro-mechanical coupling between pore pressure and exogeneous fluctuations (piezometric variations in adjacent aquifers, barometric fluctuations, and earth tides).

In addition, some aspects of the procedures developed in this work could be modified in future:

• Thus, we developed a sequential procedure to estimate first "Ss" then "Φ". First, we estimated "Ss" from wavelet analysis of semi-diurnal earth tide effects, and secondly, we estimated "Φ" from spectral analysis of barometric gain at diurnal frequency (given the previously estimated "Ss"). In reality, these 2 steps need not be implemented sequentially: they can be viewed as a system of 2 equations with 2 unknowns (Ss, Φ) to be determined both at once.

•

To estimate the global "Ss" and (also) the evolutionary "Ss(t)" during excavation, we used a fixed value, |Δε| = 2×10 -8 m 3 /m 3 , for the amplitude of the volumetric strain related to the M2 semi-diurnal earth tide (based on Melchior 1978, Bredhoeft 1967, andothers). This is because the site specific signal of volumetric strain Δε(t) was not available, so we had to estimate its amplitude indirectly from data surveys. Nevertheless, our method can be extended to the case where the strain signal Δε(t) is available. The new procedure would be to apply semi-diurnal wavelet analysis, then envelope analysis, to the Δε(t) signal, in a way similar to that already implemented for relative pore pressure p REL (t).

Another extension worth considering concerns the simplified hydromechanical "models" which we have used so far for interpreting the pressure signals. These models need to be enhanced:

• Thus, it was noted that "Ss" represents an equivalent elastic coefficient for the water filled poro-elastic rock. However, this concept rests on a simplified view of poro-elastic deformations, usually applied to confined aquifers. Its interpretation for the hydrogeologic configuration and boundary conditions of the claystone at Mont Terri is not straightforward and requires further analyses. • Also, a more "physical" hydromechanical approach based on Biot's poro-elastic theory [START_REF] Biot | General Solutions of the Equations of Elasticity and Consolidation for a Porous Material[END_REF]) could be developed, as in the coupled T-H-M models developed in [START_REF] Ababou | Thermo-Hydro-Mechanical simulation of a 3D fractured porous rock: preliminary study of coupled matrix-fracture hydraulics[END_REF] and Canamon (2009). Such models should be upgraded with explicit coupling of pore pressure to external "forcing" functions (air pressure, earth tides, etc.), in order to better interpret the effects of these signals.

•

Finally, the hydro-mechanical model could be constructed to predict space-time propagations explicitly. Thus, [START_REF] Mallet | Multidimensional modeling of a poro-elastic medium: identification of the behavior of a clay formation from a stochastic representation of piezometric fluctuations[END_REF][START_REF] Mallet | Etude des fluctuations piézométriques à l'intérieur d'un milieu poreux saturé, et détermination des paramètres caractéristiques du milieu poreux. Note Technique IRSN/DEI[END_REF]) developed a simple model to illustrate pore pressure fluctuations and propagation under boundary forcing in a compressible hydrogeologic medium.

In summary, the simplified hydro-mechanical models used in this study can be replaced, in future, with models based on PDE systems. The statistical "interpretations" of pressure signals presented in this text will then appear as the solution of PDE inverse problems (in a statistical setting). This point of view sheds a new light on the coefficients (Ss, Φ) obtained in the present study. They can be viewed as the solution of an inverse problem based on spectral gain (p REL /p ATM ) and on wavelet amplitude ratio (|Δε|/|Δp REL |). For this reason, it is thought that the statistical methods used here can be adapted for solving an inverse problem involving a hydro-mechanical PDE model forced by pressure signals and other measurements (humidity, temperature, etc.).

Acknowledgments

This chapter presents data analyses developed within Mont Terri's LP 14 experiment (Long term Pressure monitoring project, phase 14), as part of the international research programs of the MONT TERRI CONSORTIUM on the properties of claystone as a potential geologic repository for radioactive waste. The authors acknowledge the support of institutions participating in the LP 14 project at Mont Terri (ANDRA, NAGRA, IRSN, NWMO, SWISSTOPO). In addition, the authors gratefully acknowledge the help and advice provided by the following persons: Paul The cross-spectrum is complex-valued (even if the signals are realvalued).

Fig. 1 .

 1 Fig. 1. (a) Above. Location of the Mont Terri Underground Rock Laboratory in Switzerland, with its main experimental galleries. (b) Middle: Geologic cross-section of the folded Jura showing the Mont Terri anticline and the 4 km long motorway tunnel giving access to the Opalinus Clay rock (Freivogel M. 2001, Bâle University, after Schaeren G. & Norbert J. 1989: Tunnels du Mont Terri et du Mont Russelin. La traversée des "roches à risques": marnes et marnes à anhydrite. Soc. Suisse Ing. Arch., Doc. SIA D 037, 19-24). (c) Below: Schematic 3D perspective view of some of the Mont Terri galleries and the motorway tunnel.

Fig. 2 .

 2 Fig.2. Flowchart of Pre-Processing tasks and methods (programmed and implemented as a MATLAB TOOLBOX). The pre-processing tasks contain some iterations and recursions, particularly when computing statistics to detect outliers or to fill data gaps. Only some of the iterative procedures are explicitly shown here. The symbols (Mx*, Cxx*, σx*) and (mox, Coxx, σox) represent estimates of the mean, the covariance function, and the standard deviation of the signal at various stage of pre-processing. The symbol Δto is the constant time step chosen by the user for time step homogenisation. At the end of the pre-processing stage, the processed (transformed) signal is ready for statistical analyses. The processed signal may be shorter than the initial "raw" signal.

Fig. 4 .

 4 Fig. 4. Simplified flowchart presenting some of the signal analysis methods used in this work. These methods were programmed in a MATLAB TOOLBOX for SIGNAL ANALYSES.

Fig. 5 .

 5 Fig. 5. Schematic illustration of some aspects of signal sampling in discrete time within a limited observation window (which may be due to truncation of a longer time series).: t Δ time step (seconds); max t M t = Δ : size of observation window. min T : smallest period that can be identified (Shannon). max T : largest period that can be identified. max

Fig. 6 .Fig. 7 .

 67 Fig. 6. Synthetic test with 4 harmonics: (a) the RMS phase spectrum function [phase difference between X 1 (t) and X 2 (t), in radians]; (b) the RMS spectral gain function [gain between X 1 (t) and X 2 (t), in units of √X 2 /√X 1 ].

  Note: in graphs (b, c, d), each component Ci(t) is shown on a shifted ordinate axis positioned at y = i (i=0,1,2…); this allows for a better visualization of the different components on a single graph; apart from this shift, the ordinate scales are not distorted (same for all components Ci(t)).

Fig. 8 .

 8 Fig. 8. The signal X(t) and its dyadic components Ci(t) extracted from three different multiresolution bases (Haar, Daubechies 4 and Daubechies 20). (a): The synthetic non stationary signal X(t) to be analysed. (b): The decomposition of X(t) into components Ci(t) using Haar (Daubechies 2) mother wavelet. (c): The decomposition of X(t) into components Ci(t) using Daubechies 4 mother wavelet. (d): The decomposition of X(t) into components Ci(t) using Daubechies 20 mother wavelet.

Figure 9 (Fig. 9 .

 99 Figure9 (a,b) shows a construction of the successive approximations and/or residuals of signal X(t) for the lowest degree wavelet (the Haar wavelet).

Fig. 10 .

 10 Fig. 10. Signal of the raw relative pore pressure PP1(t) (kPa) at Mont Terri (from 02/08/2002 to 04/09/2002). Time span: 15 month. Time step of raw pressure signal: Δt = 30 min (and k=1).

Fig. 11 .

 11 Fig. 11. Reduced spectrum of ΔPP1(t), the first difference of relative pore pressure (kPa), at Mont Terri. Time span: ~1 month (from 02/08/2002 to 04/09/2002). Time step of raw pressure signal: Δt = 30 min (and k=1).

Fig. 12 .

 12 Fig. 12. Multiresolution wavelet analysis : several dyadic time scale components are shown for the raw relative pore pressure PP1(t) at Mont Terri. Duration: 02/08/2002 to 22/08/2002. Time step: Δt = 30 min (and sampling step: k=1).

Fig. 13 .

 13 Fig. 13. Multiresolution wavelet analysis : dyadic time scale component 8h (near 12h) for the raw signal of relative pore pressure PP1(t) at Mont Terri. Duration: 02/08/2002 to 22/08/2002. Time step: Δt = 30 min (k=1). Amplitude of fluctuations: Δh = 0.058 m.

Fig. 14 .

 14 Fig. 14. Reduced spectral gain vs. normalized frequency f×Δt (the maximum value of f×Δt is 0.5, in agreement with Shannon's sampling theorem). This spectral gain function was obtained by cross-spectral analysis of the raw signals of atmospheric pressure Patm(t) [kPa] and relative pore pressure PP1(t)-Patm(t) [kPa]. Signal duration: about 1 month (from 02/08/2002 to 04/09/2002). Time step: 30 mn. The gain at diurnal frequency (1/24h) is indicated by the arrow; it is used to estimate barometric efficiency.

Fig. 15 .Fig. 16 .

 1516 Fig. 15. Signal of the relative pore pressure PP1(t) (kPa) at Mont Terri Time span: 15 month (29/01/2004 to 12/04/2005 ). Time step of pre-processed signal: Δt = 30 min, and sampling step k=1 (i.e. no under-sampling).

  Fig. 17. Multi-resolution wavelet scale-time plot. Several dyadic components are shown here for the signal of relative pore pressure PP1(t) or p REL (t) (pre-processed) at the Mont Terri site. Signal duration: from 29/01/2004 to 12/04/2005; time step Δt = 30 min; no subsampling (k=1).

Fig. 18 .

 18 Fig. 18. Multi-resolution wavelet analysis: graph of the 8h dyadic scale component (near 12h) of the pre-processed relative pore pressure, PP1(t) or p REL (t). Record length ∼15 months (from 29/01/2004 to 12/04/2005); Δt = 30 min, k=1 (no sub-sampling). The amplitude of semi-diurnal pressure fluctuations is Δh ≈ 0.06 m equivalent water height. The methods used here are again the same as those applied in the first example of Section 4.1 where we analysed shorter, unprocessed pressure signals. We obtain here successively the specific storativity "Ss" and the effective elastic porosity "Φ". The final results are shown below for two different methods of quantifying the amplitude of semi-diurnal wavelet pressure fluctuations: |Δh| (norm 1) ≈ 0.06m; Ss (norm 1) ≈ 3.3×10 -7 m -1Φ (norm 1) ≈ 1.0% |Δh| (norm 2) ≈ 0.08m ; Ss (norm 2) ≈ 2.6×10 -7 m -1 Φ (norm 2) ≈ 0.9%

  the evolution of pore pressure and of elastic properties (specific storativity Ss; elastic porosity) 5.1 Excavation of a gallery: observed effects on the evolution of pore pressure (synexcavation phase, Mont Terri gallery Ga98)

Fig. 19 .

 19 Fig.19. Comparison of pore pressure and excavation front versus time over a 5 month period (syn-excavation). The green curve "PP1(t)" represents the time evolution of absolute pore pressure p(t) in sensor PP1. The blue curve with circles represents the distance D(t) between measurement chamber PP1 and the excavation front of gallery Ga98. The arrow at left indicates the minimal distance achieved (4 m), after which excavation was stopped for about a month (February/March 1998), before restarting around the 2 nd week of April 1998. The sharp rise of pressure PP1(t) from 2250 to 2800 kPa (Δp ≈ +5.5 bars) indicates the "fast" part of the syn-excavation phase, which lasted only 11 days (or roughly two weeks).

Fig. 20 .

 20 Fig. 20. Semi-diurnal wavelet component C 12H (t) of relative pore pressure in PP1 (blue curve) and its envelope (red curve) during the "syn-excavation" phase of gallery Ga98 at Mont Terri (5 months, from 17 Nov 1997 to 23 April 1998, with Δt=30 mn time steps).

Fig. 21 .

 21 Fig. 21. Temporal zoom of the previous figure, shown over a time window of approximately 1 month (instead of 5 months in the previous figure). The (blue) curve designated as "PP1 (12h)" is the wavelet component signal C 12H (t) extracted from PP1 relative pressure. The envelope is the (red) curve designated as "R(t)". The ∼ 11 day period of fast pressure change occurs between the dates of 7th and 19th of February 1998.

Fig. 22 .

 22 Fig. 22. Evolution of the decimal logarithm of Ss(t). The time window here is about 1.5 month or so, and it includes the 11 day period of fast change (double arrow). Here, Ss(t) is calculated from an empirical envelope (estimated with the lethod of maxima rather than the Hibert transform) of the semi-diurnal dyadic component of relative pore pressure PP1 during excavation (gallery Ga98, excavation from 17/11/1997 to 23/04/1998, time step Δt=30min).

Fig. 23 .

 23 Fig. 23. Evolution of the decimal logarithm of Ss(t), on the same time window as the previous figure (1.5 month). Here, the red curve shows Ss(t) calculated from the theoretical envelope (Hilbert Transform, Cramer-Leadbetter envelope), applied to the semi-diurnal dyadic component of relative pore pressure PP1 during excavation of gallery Ga98. The distance D(t) of the excavation front from the PP1 pressure sensor is shown as the black curve with circles "o".

Fig. 25 .

 25 Fig. 25. A broader view of the evolution of log Ss(t) (red curve, computed with the Hilbert Transform envelope). Here the time window is the entire "syn-excavation" phase of 5 months (including, again, the 11 day period of very fast changes). The black curve with circles "o" represents the distance D(t) of the excavation front from the PP1 pressure sensor. 5.5 Calculation of elastic porosity Φ(t) from barometric effect (diurnal gain p REL /p ATM ) and from the previously estimated Ss(t)

  Pore pressure (or "absolute" pore pressure), [Pa]. (*) ATM p Atmospheric pressure (or air pressure), [Pa]. (*) function of X(t) and Y(t), versus time-lag; Rxy is the normalized version of Cxy: [units of Rxy] Reduced cross-spectrum of signals X(t) and Y(t), vs. frequency "f". It is equal to the cross-spectrum

  ) in terms of continuous lag time (τ). The two parameters involved in this analytical example are the signal's correlation time (τ COR ) and the size of the observation window (t MAX ). As can be seen, the relative RMS error ε RMS is proportional to

	τ	COR MAX t	at zero lag (τ = 0), and
	ε RMS goes to 100% as

MAX t

τ →

  Bossart (SWISSTOPO); Martin Cruchaudet (ANDRA); David Jaeggi (SWISSTOPO); Herwig Mueller (NAGRA); Eric Sykes (NWMO). The authors wish also to thank Alain Mangin (CNRS-UPS Toulouse, Laboratoire Souterrain de Moulis) who initiated some of the statistical tools and methods used in this project, and others who provided scientific advice on various issues (cracks, geophysics, signals, hydrogeology, porous media): Andrea Möri (SWISSTOPO); Jacques Delay (ANDRA, LSM/HM Bure); Israel Cañamón Valera (Universidad Politecnica de Madrid); Moumtaz Razak (Université de Poitiers, HYDRASA); and Michel Quintard (IMFT, Toulouse). Reduced spectral gain vs. frequency, obtained from cross-spectrum.The complex-valued reduced gain is defined as:

	MATHEMATICAL SYMBOLS: LATIN
	Am(t)		Wavelet approximation of order "m" (multi-resolution dyadic scale
					"m")
	, B BE B , E	Barometric efficiency, [dimensionless] or [Pa/Pa] (relative pressure / atmospheric pressure effects)
	C	12 H	( ) t	Wavelet component of scale 12h (semi-diurnal)
	( ) C t m	Wavelet component of dyadic scale "m" (time scale = 2 m t)
	C	XY	( ) τ	Cross-covariance function of X(t) and Y(t), versus time-lag [units of Cxy] = [units of X(t)] [units of Y(t)]
	D(j) or D(τ(j))	Tuckey-2 filter function of time-lag (used for estimation of Fourier
					spectra)
	Dm(t)		Wavelet "detail" at dyadic scale "m" (see also: wavelet component C m (t))
	8. Appendix: Symbols and abbreviations ( ) ( ) ( ) XY XY XX g f s f s f =	.
	ABBREVIATIONS, ACRONYMS
	ANDRA	Agence Nationale pour la gestion des Déchets RAdioactifs (France)
	BPP			Borehole instrumented for Pore Pressure measurements (Mont Terri site)
	EdZ			Excavation disturbed Zone (reversibly disturbed zone around an excavation in rock, elastic hydro-mechanical disturbance)
	EDZ		Excavation Damaged Zone (irreversibly damaged, fissured and fractured zone around an excavation in rock)
	IRSN		Institut de Radioprotection et de Sûreté Nucléaire (France)
	LP14		Long term Pressure monitoring project -phase 14 (Mont Terri Project)
	NAGRA	National Cooperative for the Disposal of Radioactive Waste (Switzerland)
	NWMO	Nuclear Waste Management Organization (Canada)
	PP			Pore Pressure, i.e., absolute pore water pressure (Mont Terri site)
	SWISSTOPO	Swiss Federal Office of Topography (and Mont Terri Project)
	T-H-M or THM	Thermo-Hydro-Mechanical (coupled models)
	URL		Underground Rock Laboratory (or: Underground Research Laboratory)

f

Frequency "f" in

[Hz]

or [Number of cycles / second]; "f" is related to angular frequency "ω" by: f = ω /2π.

( )

XY g f H Total hydraulic head [m]

, related to relative pressure by:

More generally, for comparison purposes, we have also analysed pre-post-excavation phases; however the results will only be briefly mentioned here for lack of space (see Fatmi

for details).

The reader interested mostly in the applications, and not in the theoretical signal analysis methods, may jump at first reading to Sections 4-5, where the methods are implemented for Mont Terri pressure data.

These MATLAB TOOLBOXES were developed mainly by H. Fatmi, supervised by R. Ababou, with the aid of the other co-authors of this chapter and of other researchers and collaborators, including: A. Mangin, senior researcher at the Moulis CNRS laboratory (who initiated several of the signal analysis procedures); doctoral students Y. Wang and K. Alastal (who tested spectral and wavelet applications at the Institute of Fluid Mechanics); and intern students A. Mallet, A. Moulia and Ch. Joly who participated in various validation tests using synthetic signals.

The auto-covariance functions Cxx(τ) and Cyy(τ) are defined in a similar way; for instance, Cxx(τ) can be obtained simply by replacing XY by XX in the cross-covariance function Cxy(τ). The auto-correlation function Rxx(τ) is simply Cxx(τ)/Cxx(0) where Cxx(0) = Var(X(t)) = σX 2 .

Note: Shannon's sampling theorem concerns in fact only the maximum frequency; but nevertheless, there is also a minimum frequency theorem, as shown below.

Note: this example was developed recently, based on similar tests developed earlier by the authors and collaborators (see acknowledgments).

In fact, our current implementation of this general idea to the Mont Terri site is further simplified, using constant tidal amplitude because of the lack of direct earth dilatation record at the site.

Cramer-Leadbetter theory in a nutshell: see Box N°1.

The (auto)-spectrum is always real-valued for a real-valued process X(t).

S S

Specific elastic storativity, [m -1 ] or [m 3 /m 3 per meter of head δH] (*) Pressures "p" are sometimes denoted with a capital "P" within the text.

MATHEMATICAL SYMBOLS: GREEK

, TIDAL ε ε Δ Δ Tidal dilatation, volumetric strain [m 3 /m 3 ] (here, this notation is used for the amplitude of the tidal strain).