Elizaveta Saifutdinova

Marco Congedo
email: marco.congedo@gipsa-lab.fr

Daniela Dudysova
email: daniela.dudysova@nudz.cz

Lenka Lhotska
email: lenka.lhotska@cvut.cz

Jana Koprivova
email: jana.koprivova@nudz.cz

Vaclav Gerla
email: vaclav.gerla@cvut.cz

An Unsupervised Multichannel Artifact Detection Method for Sleep EEG Based on Riemannian Geometry

Keywords: sleep EEG, artifact detection, Riemannian geometry

In biomedical signal processing, we often face the problem of artifacts that distort the original signals. This concerns also sleep recordings, such as EEG. Artifacts may severely affect or even make impossible visual inspection, as well as automatic processing. Many proposed methods concentrate on certain artifact types. Therefore, artifact-free data are often obtained after sequential application of different methods. Moreover, single-channel approaches must be applied to all channels alternately. The aim of this study is to develop a multichannel artifact detection method for multichannel sleep EEG capable of rejecting different artifact types at once. The inspiration for the study is gained from recent advances in the field of Riemannian geometry. The method we propose is tested on real datasets. The performance of the proposed method is measured by comparing detection results with the expert labeling as a reference and evaluated against a simpler method based on Riemannian geometry that has previously been proposed, as well as against the state-of-the-art method FASTER. The obtained results prove the effectiveness of the proposed method.

Introduction

Human electroencephalography (EEG) refers to recordings of brain activity obtained using electrodes placed on the surface of the head. EEG investigation is one of the main research methods in sleep medicine. It aims at recording electrical activity caused by internal processes during pre-sleep wakefulness and sleep. The number of electrodes and their locations are set according to the areas of interest. Modern sleep research often utilizes at least 19 channels distributed evenly over the head. Some recent studies have applied up to 256-channel high-density EEG [START_REF] Siclari | Dreaming in NREM Sleep: A High-Density EEG Study of Slow Waves and Spindles[END_REF][START_REF] Ratti | A high-density polysomnographic picture of disorders of arousal[END_REF]. Manual data inspection is a demanding and tedious task, especially in long-term recordings such as a whole-night EEG. This makes automatic processing methods highly needed in sleep medicine research.

Events not related to brain activity may noticeably affect automatic analysis or make it even impossible. Such events are called artifacts. Typical causes of artifacts are technical issues (detached electrodes, sweating, electrical noise) and physiological events (muscle contractions, movements, eye rolling) of the body. For instance, artifacts caused by blinking are characterized by slow-frequency/high-amplitude waves that can be seen mainly in frontal electrodes. The typical waveform corresponding to such artifacts is shown in Figure 1a. Blinking occurs only during the awake state. Artifacts associated with movements may arise during both being awake and asleep, the latter divided into rapid eye movement (REM), non-REM 1, 2, and 3 sleep stages [START_REF] Iber | American Academy of Sleep Medicine. The AASM Manual for the Scoring of Sleep and Associated Events: Rules, Terminology and Technical Specifications[END_REF]. Such artifacts are characterized by an unusually high amplitude and unusual frequency in several electrodes, as is shown in Figure 1b. Electrode artifacts may arise during the whole recording, and they may have different characteristics. Voltage jumps could be seen in all channels and are illustrated in Figure 1c. A loss of contact of electrodes may produce noise in a channel or high amplitude anomalies, as shown in Figure 1d. Manual artifact detection is a time-consuming task requiring specific skills. Manual detection is more difficult in high density EEG recordings (128 or 256 channels), and it might lead to a higher number of undetected artifacts on one side and to false detection on the other side. There is a wide choice of automated artifact detection methods available in the literature [START_REF] Islam | Methods for artifact detection and removal from scalp EEG: A review[END_REF]. They include approaches based on signal decomposition such as Fourier transform [START_REF] Szentkirályi | Performance of an automated algorithm to process artefacts for quantitative EEG analysis during a simultaneous driving simulator performance task[END_REF], wavelet transform [START_REF] Taherisadr | Single Channel EEG Artifact Identification Using Two-Dimensional Multi-Resolution Analysis[END_REF][START_REF] Anastasiadou | Unsupervised detection and removal of muscle artifacts from scalp EEG recordings using canonical correlation analysis, wavelets and random forests[END_REF], or empirical mode decomposition [START_REF] Bono | Hybrid wavelet and EMD/ICA approach for artifact suppression in pervasive EEG[END_REF][START_REF] Chen | The Use of Multivariate EMD and CCA for Denoising Muscle Artifacts From Few-Channel EEG Recordings[END_REF]. A popular method for eye blinking artifact detection is regression [START_REF] Klados | ICA: A hybrid methodology combining Blind Source Separation and regression techniques for the rejection of ocular artifacts[END_REF][START_REF] Borowicz | Using a multichannel Wiener filter to remove eye-blink artifacts from EEG data[END_REF]. Blind source separation (BSS) methodologies like independent component analysis (ICA) [START_REF] Islam | Methods for artifact detection and removal from scalp EEG: A review[END_REF][START_REF] Szentkirályi | Performance of an automated algorithm to process artefacts for quantitative EEG analysis during a simultaneous driving simulator performance task[END_REF][START_REF] Mayeli | Real-time EEG artifact correction during fMRI using ICA[END_REF] and independent vector analysis [START_REF] Chen | Independent Vector Analysis Applied to Remove Muscle Artifacts in EEG Data[END_REF] are often used for ocular, muscular, and cardiac artifact rejection. Once the decomposition has been achieved, components corresponding to artifact activity may be rejected manually or using automatic methods [START_REF] Chang | Evaluation of Artifact Subspace Reconstruction for Automatic EEG Artifact Removal[END_REF]. A combination of BSS and other signal decomposition techniques may bring benefits for muscular artifact removal [START_REF] Chen | A Novel EEMD-CCA Approach to Removing Muscle Artifacts for Pervasive EEG[END_REF][START_REF] Chen | A Preliminary Study of Muscular Artifact Cancellation in Single-Channel EEG[END_REF]. Many recent studies have proposed methods based on classification techniques such as random forest classifiers [START_REF] Anastasiadou | Unsupervised detection and removal of muscle artifacts from scalp EEG recordings using canonical correlation analysis, wavelets and random forests[END_REF], artificial neural networks [START_REF] Nguyen | EOG artifact removal using a wavelet neural network[END_REF], or support vector machine classifiers [START_REF] Levitt | Automated detection of electroencephalography artifacts in human, rodent and canine subjects using machine learning[END_REF][START_REF] Dhindsa | Filter-Bank Artifact Rejection: High performance real-time single-channel artifact detection for EEG[END_REF], where training examples determine the types of artifacts to be detected. Another direction of modern artifact detection is the optimization of existing approaches [START_REF] Acharyya | Low-complexity hardware design methodology for reliable and automated removal of ocular and muscular artifact from EEG[END_REF][START_REF] Gerla | Automatic identification of artifacts and unwanted physiologic signals in EEG and EOG during wakefulness[END_REF], adjusting them to the increasing volume of recorded data. Research in the field of sleep data processing mainly focuses on the detection and elimination of the most expected artifacts during sleep [START_REF] Wallant | Automatic artifacts and arousals detection in whole-night sleep EEG recordings[END_REF] through single-channel approaches [START_REF] Malafeev | Automatic artefact detection in single-channel sleep EEG recordings[END_REF]. These methods concentrate only on certain artifact types. Therefore, artifact-free data are obtained after sequential application of different methods, and single-channel approaches must be applied to all channels alternately.

The aim of the study is to develop an artifact detection method in multichannel sleep EEG capable of rejecting different artifact types at once. The inspiration for the study is gained from the idea of the Riemannian potato [START_REF] Barachant | The Riemannian Potato: An Automatic and Adaptive Artifact Detection Method for Online Experiments Using Riemannian Geometry[END_REF], an automatic and adaptive artifact detection method proposed for online brain-computer interface (BCI) experiments. It considers spatial patterns of awake EEG data and exploits the fact that clean data tend to form a cluster in the manifold of symmetric positive matrices (the "potato"), the space on which the EEG data is mapped, while artifacts tend to be distributed outside the cluster. In our study, we assume that the number of clusters for clean data may be superior to one, and we introduce an unsupervised and adaptive method for constructing the clusters. It allows us to achieve better performance for complex data such as sleep EEG, where we cannot assume a unique cluster because the spatial patterns vary widely across sleep stages and because sleep data within the same stage may be captured more accurately by several clusters. Different sleep stages have their own characteristics. In the present study, we investigate the proposed method and test it on real sleep EEG data. Performance is evaluated by comparing the obtained detection to a provided expert scoring. Results are examined in contrast to the outcome of the (one-cluster) Riemannian potato artifact detector [START_REF] Barachant | The Riemannian Potato: An Automatic and Adaptive Artifact Detection Method for Online Experiments Using Riemannian Geometry[END_REF] and the fully-automated statistical thresholding for EEG artifact rejection method (FASTER) [START_REF] Nolan | Fully Automated Statistical Thresholding for EEG artifact Rejection[END_REF], which performs artifact detection by analyzing statistical properties in separate channels.

Data

Dreams Dataset

The open-source Dreams dataset [START_REF] Devuyst | Automatic Processing of EEG-EOG-EMG Artifacts in Sleep Stage Classification[END_REF] contains 15-min recordings of sleep data, 20 recordings in total. Channels originally referenced to A1 were re-referenced to Cz. All available EEG channels O1-Cz, O2-Cz, Fp1-Cz, and Fp2-Cz were used in the study. The sampling frequencies were 50 Hz, 100 Hz, and 200 Hz. We excluded from analysis the recordings sampled at 50 Hz due to the fact that with such a sampling rate, the frequency band-pass is too restricted for sleep data analysis [START_REF] Iber | American Academy of Sleep Medicine. The AASM Manual for the Scoring of Sleep and Associated Events: Rules, Terminology and Technical Specifications[END_REF]. Furthermore, we excluded recording Numbers 3, 6, 10, and 16 because the proportion of artifacts exceeded 40%. Sleep stage scoring and artifact labeling have been performed by a trained clinician and provided, as well. Each recording contains a mix of sleep stages. The distribution of the sleep stages in the dataset is: wakefulness 33%, REM 25%, non-REM 1 18%, non-REM 2 23%, and other 1%. Information on the artifact rate is shown in Table 1.

InSleep Datasets

The InSleep datasets were obtained from the National Institute of Mental Health, Czech Republic. The study was conducted in accordance with the Declaration of Helsinki, and the protocol was approved by the Ethics Committee of the National Institute of Mental Health (Project Number 6/15). The whole-night EEG as a part of the polysomnography (PSG) was recorded from subjects suffering from insomnia. Each recording was obtained using 19 channels placed according to the 10-20 system: Pz, Cz, Fz, T6, T5, T4, T3, F8, F7, O2, O1, P4, P3, C4, C3, F4, F3, Fp2, Fp1. Channels were referenced to an average of Cz, Fz, and Pz. The sampling frequency was 250 Hz. Sleep stages were scored by a trained clinician. Each recording contained 10-20 min of data corresponding to a single sleep stage (REM, non-REM 2, or non-REM 3) or awake activity. The non-REM 1 sleep stage was not presented in the dataset due to its short duration in sleep recordings. Artifacts were labeled by an expert. Visual inspection of EEG was performed separately for each recording in the EEG Lab software using function "reject continuous data by eye" [START_REF] Delorme | An open source toolbox for analysis of single-trial EEG dynamics including independent component analysis[END_REF]. The expert was instructed to label all visible artifacts. The total number of obtained recordings was 44. Details of the artifact rate are provided in Table 1.

Materials and Methods

Theoretical Background

This section describes the theoretical aspects of working with multichannel data. We denote by x t ∈ R n the signal vector at n channels and at time point t. X = [x s , ..., x s+m-1 ] ∈ R n×m denotes an epoch starting at sample s and lasting m samples. For X, we estimate the spatial covariance matrix using the usual sample covariance matrix (SCM) estimator:

C = 1 m -1 XX T ∈ R n×n . (1)
The matrix C is symmetric positive-definite (SPD). Often, SPD matrices are considered in a Euclidean space with the associated Frobenius norm C F . However, the native space of SPD matrices is not Euclidean, but a Riemannian space [START_REF] Congedo | Riemannian geometry for EEG-based brain-computer interfaces; a primer and a review[END_REF]. Riemannian geometry starts by defining an inner product (metric) at each tangent space to the manifold, which varies smoothly from point to point. Applying the affine-invariant (Fisher) metric [START_REF] Congedo | Riemannian geometry for EEG-based brain-computer interfaces; a primer and a review[END_REF], the distance δ between two points C i and C j on the manifold is given by:

δ(C i , C j ) = n ∑ i=1 log 2 λ i 1/2 , ( 2
)
where the λ i are the eigenvalues of C -1 i C j . A centroid (called the geometric mean) of an SPD matrix set can be defined using the above distance [START_REF] Congedo | A fixed-point algorithm for estimating power means of positive definite matrices[END_REF]. The geometric mean M of a covariance matrix set C i is calculated using an iterative algorithm. First, M is initialized by the arithmetic mean of C i . Then, the following iteration is repeated:

M = M 1/2 exp[1/k ∑ k ln(M -1/2 C k M -1/2 )]M 1/2 (3)
until convergence, which is obtained when:

∑ k ln(M -1/2 C k M -1/2 ) F < ε. ( 4
)
The distance of a set of covariance matrices to their geometric mean does not have a symmetric distribution. To make the distribution symmetric, we may use the standardized distances to the geometric mean [START_REF] Congedo | EEG Source Analysis[END_REF]; define the geometric mean µ of the distances δ k as:

µ = exp( 1 k ∑ k ln δ k ) (5)
and the geometric standard deviation σ of the distances as:

σ = exp( 1 k ∑ k (ln δ k µ ) 2 ). (6)
Then, the standardized distance measures are:

δ k = ln δ k µ ln σ . ( 7
)

The Riemannian Potatoes' Artifact Detection

A Riemannian potatoes artifact detection method (RPs) is proposed in this study. It relies on the assumption that on the SPD manifold, contaminated epochs are mapped far away from clusters of clean data. In contrast to the simple potato [START_REF] Barachant | The Riemannian Potato: An Automatic and Adaptive Artifact Detection Method for Online Experiments Using Riemannian Geometry[END_REF], here we take into consideration the fact that several clusters may be needed to describe the distribution of clean sleep data. In the first step, RPs restore these clusters from a recorded EEG. In the second step, we obtain continuous scoring by analyzing the data in a sliding window and comparing them to the obtained clusters. Finally, artifacts are detected by applying a threshold. The method overview is depicted in Figure 2. As in [START_REF] Barachant | The Riemannian Potato: An Automatic and Adaptive Artifact Detection Method for Online Experiments Using Riemannian Geometry[END_REF], we filtered the data before further processing. However, we used a low-pass filter under 30 Hz to save frequencies significant for sleep stages [START_REF] Iber | American Academy of Sleep Medicine. The AASM Manual for the Scoring of Sleep and Associated Events: Rules, Terminology and Technical Specifications[END_REF]. The cluster construction process starts with splitting the filtered data into non-overlapping 1-s epochs. Such an epoch length is a standard for Riemannian geometry analysis [START_REF] Barachant | The Riemannian Potato: An Automatic and Adaptive Artifact Detection Method for Online Experiments Using Riemannian Geometry[END_REF]. Covariance matrices are extracted from every epoch. Then, distances between all pairs of covariance matrices are calculated by Equation (2). The average distance d i is the mean value of distances from the ith epoch to all the other epochs. Epochs with a d value greater than a threshold are removed from further analysis. We consider the mean value of d as the threshold. Then, we cluster the remaining matrices using the k-means algorithm with the distance measure defined by Equation (2). To determine the number of clusters, we use the following procedure: the number of clusters increases iteratively from one to a maximum of 10; centroids are obtained as the geometric means; the process is finished once the constructed clusters satisfy an overall normality condition, that is when the distribution of standardized distances to the centroids is normal in all clusters. Normality is tested using D'Agostino's K 2 test. The distribution of distances to the centroid within a cluster is examined separately. The overall normality p-value is obtained as a combination of the obtained p-values using Stouffer's Z-score method. The cluster construction process is stopped when the combined p-value is superior to 0.05. If the cluster construction finishes without achieving overall normality, then we choose the number of clusters achieving the greatest p-value. In the case of tied p-values, we retain the one corresponding to the smaller number of clusters. The rationale for this procedure is that when the distances of the points in the cluster from the center of mass have a Gaussian distribution, we assume that the cluster is a good representative of a region in the manifold covered by the data and that further splitting is therefore unnecessary.

In the second step, RPs analyze the n-channel data in a sliding window of a size of 1 s with a step equal to 0.1-times the sampling frequency. Assume there are k clusters obtained in the previous step and denote O j ∈ R n×n , j ∈ [1, k] their centroids. For each ith window, we obtain covariance matrix C i ∈ R n×n . All k distances σ ij between C i and every O j centroid are calculated. Then, the closest cluster j * to C i is determined by the minimum distance among the σ ij values. The geometric mean and standard deviation within the cluster were determined in the cluster construction step. Therefore, we can obtain σ * i as a standardized distance σ ij * to the centroid O j * by Equation [START_REF] Anastasiadou | Unsupervised detection and removal of muscle artifacts from scalp EEG recordings using canonical correlation analysis, wavelets and random forests[END_REF]. Due to the cluster construction, standardized distances within a cluster are distributed normally. We can consider the obtained σ * i as a z-score and translate it into a probability value p i using the Gaussian cumulative distribution function. The value p i indicates a probability of epoch i with covariance matrix C i to belong to a given cluster. Hence, the value of p i equals 1p i and indicates a probability of epoch i with covariance matrix C i to be an outlier.

We consider that the values σ * i correspond to the center of the ith epoch and that in an epoch, the scoring starts and ends with 0. Therefore, we interpolate the obtained values within the epoch. The obtained score is then smoothed with a moving average filter to make the output more robust. The local minima of these scores define segments. Finally, a segment whose score is above a threshold for at least 0.4 s is considered as an artifact. This value was established as a minimum artifact duration by the experts.

Experimental Results

Results Evaluation

Artifact detection performance was evaluated by comparing the detection obtained by the automatic artifact detection methods to a scoring provided by human experts. The comparison was performed on a sample-by-sample basis. We denote the number of true positives, false positives, true negatives, and false negatives as TP, FP, TN, and FN respectively.

Precise artifact borders are hard to determine. Human and automatic scoring may differ and yet be both considered correct, as is shown in Figure 3. To handle this, we added a fuzzy area on each artifact border labeled by an expert. Thus, residual segments with no matching labeling as Segments A-B and C-D in Figure 4 were treated as TP and TN, respectively, if they lasted for at least 10% of the artifact duration and no more than 1.5 s. Otherwise, they were counted as FP or FN, respectively. The parameters used for results evaluation were estimated based on experience and consultancy with neurological experts who described their visual inspection process. Following [START_REF] Wallant | Automatic artifacts and arousals detection in whole-night sleep EEG recordings[END_REF], we computed Cohen's kappa K and the percentage of agreement S. K is considered a more robust measure than S because it takes into account agreement occurrence by chance. Values of K greater than 0.6 are considered as good and values higher than 0.8 indicate an excellent agreement [START_REF] Landis | The Measurement of Observer Agreement for Categorical Data[END_REF]. We have:

K = P o -P r 1 -P r , (8)
where:

P o = TP + TN TP + TN + FP + FN , (9)
P r = (TP + FN) * (TP + FP) + (FP + TN) * (TN + FP) (TP + TN + FP + FN) 2 . ( 10
)
Sensitivity Se and false discovery rate (FDR) were also compared in this study.

Se = TP TP + FN , (11)
FDR = FP TP + FP . ( 12
)
The non-parametric Wilcoxon signed rank test was used for comparing metrics. This test does not require normality of the samples and is also suitable for small sample sizes.

Results

Our RPs artifact detection method was implemented in Python and applied to the Dreams and InSleep datasets. The Riemannian potato [START_REF] Barachant | The Riemannian Potato: An Automatic and Adaptive Artifact Detection Method for Online Experiments Using Riemannian Geometry[END_REF] method was chosen as the benchmark. We used the original implementation of the Riemannian potato artifact detector provided by the authors in the covariance toolbox for MATLAB (github.com/alexandrebarachant/covariancetoolbox). Visual comparison of scoring obtained by RPs and by the benchmark is provided in Figures 5 and6. Detection results were compared with scoring provided by an expert, and the three statistical metrics presented in the previous section were computed as well as percentage of data labeled as artifacts by methods. Each recording of the Dream dataset contains a mix of sleep stages and awake EEG. Some of the presented sleep stages have a short duration and do not include artifacts. Testing results grouped for such EEG types would not be representative. Therefore, metric values for the Dreams dataset were obtained for the entire recordings. Recordings in the InSleep dataset contain activity of a single stage. That allowed us to access the method performance on separate sleep stages and wakefulness. Table 2 reports the obtained results for the benchmark, and Table 3 contains values for the proposed method. Cohen's kappa for both methods on all datasets is shown in Figure 7. The proposed method, RPs, performed better than the benchmark in terms of FDR for the Dreams dataset (p-value < 0.05) and InSleep REM and non-REM 2 data (p-value < 0.05 and p-value < 0.01, respectively). With respect to sensitivity, RPs did not outperform the benchmark. Finally, RPs significantly outperformed the benchmark method on non-REM 2 (p-value < 0.01) and non-REM 3 (p-value < 0.05) in terms of agreement with the expert scoring. A similar trend is observed for the InSleep dataset, as shown in Figure 7 (p-value < 0.1). Visual inspection in the Dreams dataset revealed that the benchmark method tended to detect incorrectly as artifacts characteristic EEG patterns of non-dominant sleep stages such as high amplitude delta-waves corresponding to non-REM 3 activity in recordings with a dominant non-REM 2 sleep stage, as illustrated in Figure 8. FDR values obtained on InSleep non-REM 2 indicated a big proportion of false positive events. Both methods incorrectly detected injections of non-REM 3 sleep in non-REM 2. This activity is defined as high amplitude slow waves lasting 1-5 s. An example is provided in Figure 9a. Such events account for a small proportion of the recording for which the methods do not appear adapted. Furthermore, the non-REM 2 sleep data also contained K-complexes and sleep spindles characterized as a single high-amplitude delta wave and sigma bursts, respectively. An example is depicted in Figure 9b. A wide range of patterns of such events complicates the detection for both methods. Moreover, artifacts in stages non-REM 2 and non-REM 3 occurred rarely according to Table 1; this made the metrics values for these datasets more sensitive to false detections. As we have mentioned, RPs utilizes 1-s epochs, which is a standard window length for Riemannian geometry analysis. We have also applied RPs using 2-s and 4-s epoch lengths to all datasets to test whether longer epoch lengths are more suitable for capturing the dynamics of slow waves. The obtained results are presented in Tables 4 and5 respectively. The results showed no significant increase in K in all datasets. In fact, a significant decrease in K was observed for wakefulness and REM InSleep (p-value < 0.01). Indeed, the number of events corresponding to normal slow wave activity and incorrectly determined as artifacts in non-REM 2 and 3 decreased. However, using such modifications decreased the number of detected events in general and increased the number of false results. This is caused by an increased number of missed small artifacts. Moreover, the detected events were too long, which led to an increase of the FP value and, consequently, FDR.

time, s

Additionally, we have applied FASTER [START_REF] Nolan | Fully Automated Statistical Thresholding for EEG artifact Rejection[END_REF] to all datasets. This method rejects epochs of a single-channel EEG based on calculated thresholds. The process is repeated for each channel. Mean, variance, amplitude range, and mean gradient are used for decision making. The results are presented in Table 6. The proposed RPs method significantly outperformed FASTER in K for all datasets (p-value < 0.05 for InSleep non-REM 3, p-value < 0.01 for others). Moreover, RPs achieves significantly better Se in InSleep wakefulness (p-value < 0.01), REM, and non-REM 2 (p-value < 0.05 for both datasets) and significantly smaller FDR for Dreams, InSleep non-REM 2 and 3 (p-value < 0.01, p-value < 0.05, and p-value < 0.01, respectively). The method performed well in the detection of movement artifacts. However, often, it incorrectly detected a high amplitude delta in non-REM 2 (Figure 9a) and non-REM 3 (Figure 9c) as an artifact. Furthermore, it missed many short and small-amplitude artifacts (Figure 9d). For the estimation of time consumption, we ran FASTER and RPs on 15-min 19-channel recordings. A Windows 10 Pro computer with an Intel Core i7 2.40-GHz CPU and 16 GB of RAM was used. Testing of FASTER, excluding the ICA decomposition, which was very time consuming, was performed using MATLAB R2014b. The processing time was 5.44 ± 0.01 s per recording. Testing of RPs was performed using Python 2.7. The processing time of RPs was 43.86-112.63 s per recording depending on the number of constructed clusters. Cluster construction took 0.03 s for a single cluster and 68.72 s for the construction of six clusters. Testing on the Dreams dataset resulted in 1-4 clusters. The number of clusters in the InSleep dataset varied in the range of 1-6 for awake and REM EEG and in the range of 1-4 for non-REM 2 and 3 activity. The number of clusters did not depend on the number of channels, but only on the presence of different activities, which could be distinguished in terms of the Riemannian geometry of the data.

Discussion

The performance of RPs was compared to a previously-developed artifact detection method based on Riemannian geometry. Testing on rest-state awake data confirmed the effectiveness of both RPs and the benchmark approaches. The results obtained for sleep data proved that the proposed method was more favorable. The reason why RPs outperformed the benchmark in sleep data is that the multiple cluster construction was adapted to the variability of normal spatial patterns in sleep stages. This resulted in a reduction of false detections and improved agreement with the expert. The method was fully unsupervised and adaptive; therefore, detection was dependent only on the input data and could be completely automated. Testing results proved that on the data we considered, a 1-s epoch length was an optimal choice for RPs. Furthermore, the proposed method is scalable and may be applied for any number of channels. Nevertheless, using highly-correlated signals makes the application of the method cumbersome because the estimated covariance matrices may be badly conditioned. Therefore, in the case of a very large number of EEG channels, a dimensionality-reduction pre-processing step is recommended. To do this, methods inspired by the Riemannian geometry of SPD matrices may be used; see [START_REF] Congedo | A closed-form unsupervised geometry-aware dimensionality reduction method in the Riemannian Manifold of SPD matrices[END_REF][START_REF] Rodrigues | Dimensionality Reduction for BCI classification using Riemannian geometry[END_REF].

In comparison to FASTER, which relies on statistical properties in separate channels, the methods based on Riemannian geometry achieved better performance on all datasets. FASTER missed too many artifacts, mostly small ocular and short electrode artifacts. Nevertheless, all tested methods equally performed well on the detection of high frequency/amplitude artifacts. Such artifacts prevailed in wakefulness, which allowed FASTER to achieve better performance among all datasets therein. For sleep, FASTER overestimated artifact activity, which lad to low performance. In non-REM 2, it detected only injections of non-REM 3 activity, which noticeably affected the metric values. All methods detected injections of non-REM 3 sleep into non-REM 2 as an artifact. Even if such events were not labeled by the expert as artifacts, their elimination led to improvement of the spectral characteristics of non-REM 2 activity. As for computational efforts, FASTER significantly outperformed RPs, since it utilized simple signal statistics, whereas RPs was based on eigenvalue estimations for computing Riemannian distances. The time required by the RPs method, however, is low in absolute terms and very acceptable for any practical purpose.

This study shows once more that automatic artifact detection in sleep EEG is a challenging task and demonstrates the low efficiency of methods developed for awake EEG. Many studies addressing this problem have proposed methods for artifact detection without information about sleep stages [START_REF] Wallant | Automatic artifacts and arousals detection in whole-night sleep EEG recordings[END_REF]. However, EEG investigation pipelines in modern sleep research laboratories include analysis of sleep data with already labeled sleep stages. Some recent studies have proposed artifact detection based on analysis of 20-and 30-s epochs [START_REF] Malafeev | Automatic artefact detection in single-channel sleep EEG recordings[END_REF]. Evaluation using 20-s (30-s) epochs is required for data investigation at a macro level in sleep scoring by humans, whereas for an automatic artifact rejection, it can be performed at a micro-level, and it is preferable to consider smaller epochs. A way to make the automatic method more comparable to human scoring would be to reject whole 20-s (30-s) segments if the majority of 1-s epochs forming them were rejected. Moreover, sleep disorders increase the inter-subject variability of the EEG, as well as the number of artifacts, which makes the task of artifact detection even more demanding. For these and other situations, a fully-unsupervised and adaptive strategy such as RPs is theoretically appealing.

Limitations and Future Work

The RPs method we have proposed has several disadvantages and therefore can be further improved. First, infrequent EEG patterns might not be represented in the clusters of artifact-free activity due to the cluster construction procedure. Thus, they might be incorrectly detected as artifacts.

An improvement could be achieved by setting an automatic method for identification of expected patterns like K-complexes and post-processing analysis. Moreover, the method seems to fail in rejecting artifacts that are not widely distributed on the scalp, such as muscular or cardiac artifacts. This may be obviated by constructing smaller potatoes that include only a small number of channels and by combining the rejection of all the potatoes (Riemannian potato fields), which is currently under investigation. The method was not tested on high-density EEG. In this case, the problem of highly-correlated signals in high-density EEG must be carefully considered. Future work will also be focused on improvements of the method in order to eliminate false detections of K-complexes and sleep spindles. The application of the method is currently under consideration for semi-automatic artifact rejection at the National Institute of Mental Health of Czech Republic. At this stage, probability scoring will be provided along with data to support manual investigation. Additionally, the method will be extended to multimodal PSG recordings. Other directions of future research include application of the proposed method and analysis of the obtained clusters for automatic sleep stage identification using supervised machine learning approaches.

Conclusions

The study presents an unsupervised multichannel artifact detection method for sleep EEG based on Riemannian geometry. The presented method forms clusters of artifact-free EEG data considering their spatial patterns. Application of the proposed method to sleep recordings brings significant benefits. The method is scalable, fully unsupervised and adaptive, independent of artifact types, and the outcome is easily interpreted. In comparison to the Riemannian potato artifact detector, it demonstrates better performance on complex sleep data in terms of agreement with human scoring and reduces the number of events incorrectly identified as artifacts. The RPs toolbox is available at https://gitlab.ciirc.cvut.cz/open-source/rps.

Figure 1 .

 1 Figure 1. Expert labeling (Expert) of blinking artifacts (a), movement artifact during wakefulness (b), and electrode artifacts in non-REM 3 (c) and REM (d). Solid blue lines are the expert's artifact marks.

Figure 2 .

 2 Figure 2. The method overview.

Figure 3 .

 3 Examples of comparison of human scoring and artifact detection (a,b).

Figure 4 .

 4 Figure 4. Detection evaluation. Horizontal bold lines represent detected segments. The horizontal dotted line stands for segments labeled as artifact-free data. Vertical dashed lines are the ends of all labeled segments.

Figure 5 .

 5 Figure 5. Scoring obtained by the benchmark method (RP) and RPs (RPS) on data from the Dreams dataset. Normalized distances to four clusters obtained for this recording with RPs (CL1-4) are shown, as well. Artifacts (Expert) are denoted with the solid blue line.

Figure 6 .Figure 7 .

 67 Figure 6. Scoring obtained by the benchmark method (RP) and RPs (RPS) on data from the InSleep dataset. Normalized distances to three clusters obtained for this recording with RPs (CL1-3) are shown, as well. Artifacts (Expert) are denoted with a solid blue line.

Figure 9 .

 9 Figure 9. Expert scoring (Expert), results of the FASTER (FASTER), benchmark (RP), and the proposed (RPS) methods of different EEG patterns in delta waves (a), k-complex in non-REM 2 (b), delta waves in non-REM 3 (c), and ocular artifact in REM (d).

Table 1 .

 1 Dataset details.

Dataset Sleep Stage Number of Subjects Length (minutes) Artifact Rate (%)

	Dreams		12	15	9.38 ± 9.84
		Wakefulness	14	15.38 ± 2.26	18 ± 9.9
	InSleep	REM Non-REM 2	12 10	13.5 ± 1.88 15.3 ± 1.77	7.43 ± 5.82 4.07 ± 5.33
		Non-REM 3	8	14.75 ± 2.76	1.97 ± 1.43

Table 2 .

 2 Performance of the benchmark artifact detection method.

	Dataset Sleep Stage	K	Se	FDR	Artifact Rate (%)
	Dreams		0.58 ± 0.2	0.79 ± 0.19	0.4 ± 0.21	12.23 ± 13.71
		Wakefulness 0.84 ± 0.09 0.91 ± 0.08 0.16 ± 0.07	15.2 ± 8.24
	InSleep	REM Non-REM 2 0.40 ± 0.24 0.77 ± 0.21 0.62 ± 0.41 0.75 ± 0.24 0.88 ± 0.12 0.28 ± 0.24	8.3 ± 7.14 9.12 ± 13.53
		Non-REM 3 0.54 ± 0.17 0.62 ± 0.29 0.37 ± 0.21	1.58 ± 1.52

Table 3 .

 3 Performance of RPs artifact detection.

	Dataset Sleep Stage	K	Se	FDR	Artifact Rate (%)
	Dreams		0.74 ± 0.17 0.78 ± 0.18 0.24 ± 0.22	7.49 ± 7.65
		Wakefulness 0.85 ± 0.09	0.9 ± 0.08	0.12 ± 0.08	16.6 ± 7.62
	InSleep	REM Non-REM 2 0.50 ± 0.37 0.84 ± 0.17 0.55 ± 0.38 0.82 ± 0.11 0.86 ± 0.12 0.18 ± 0.13	9.47 ± 8.22 9.06 ± 5.81
		Non-REM 3 0.68 ± 0.22 0.68 ± 0.26 0.22 ± 0.22	1.42 ± 1.5

Table 4 .

 4 Performance of RPs artifact detection with a 2-s epoch length. REM 2 0.53 ± 0.41 0.79 ± 0.21 0.49 ± 0.44 9.75 ± 7.78 Non-REM 3 0.67 ± 0.26 0.67 ± 0.26 0.26 ± 0.27 2.74 ± 2.08

	Dataset Sleep Stage	K	Se	FDR	Artifact Rate (%)
	Dreams		0.71 ± 0.17	0.85 ± 0.1	0.31 ± 0.17	12.17 ± 13.29
		Wakefulness 0.77 ± 0.13	0.85 ± 0.1	0.2 ± 0.14	18.75 ± 8.17
	InSleep	REM Non-	0.72 ± 0.14	0.8 ± 0.19	0.25 ± 0.19	8.99 ± 5.31

Table 5 .

 5 Performance of RPs artifact detection with a 4-s epoch length.

	Dataset Sleep Stage	K	Se	FDR	Artifact Rate (%)
	Dreams		0.52 ± 0.16	0.8 ± 0.14	0.54 ± 0.16	15.12 ± 15.27
		Wakefulness 0.68 ± 0.11	0.84 ± 0.1	0.32 ± 0.14	20.75 ± 8.62
	InSleep	REM Non-REM 2 0.46 ± 0.38 0.77 ± 0.19 0.57 ± 0.38 0.53 ± 0.16 0.66 ± 0.2 0.46 ± 0.2	9.59 ± 7.92 9.52 ± 7.74
		Non-REM 3 0.55 ± 0.23	0.62 ± 0.3	0.46 ± 0.24	2.87 ± 2.18

Table 6 .

 6 Performance of the FASTER artifact detection method.

	Dataset Sleep Stage	K	Se	FDR	Artifact Rate (%)
	Dreams		0.31 ± 0.22 0.81 ± 0.28 0.66 ± 0.26	24.9 ± 11.55
		Wakefulness 0.58 ± 0.21 0.52 ± 0.24 0.09 ± 0.14	6.45 ± 6.14
	InSleep	REM Non-REM 2 0.28 ± 0.25 0.43 ± 0.36 0.72 ± 0.25 0.46 ± 0.21 0.57 ± 0.24 0.45 ± 0.27	8.64 ± 6.3 3.46 ± 2.12
		Non-REM 3 0.41 ± 0.21 0.66 ± 0.34 0.64 ± 0.16	3.29 ± 1.57

Author Contributions: E.S. was co-responsible for the development of the proposed methods, program coding, paper writing, and revision. M.C. was co-responsible for development of the proposed methods, participated in the study conception and design, and provided a review. D.D. provided the data collection and labeling. L.L. provided funding acquisition and participated in drafting and critical revision of the manuscript. J.K. participated in drafting and critical revision of the manuscript. V.G. conducted research for the related works and helped with data analysis.

Conflicts of Interest:

The authors declare no conflict of interest.

Abbreviations

The following abbreviations are used in this manuscript:


Eva Moreno
email: eva.moreno@mnhn.fr

Catherine Homberg

Johann Schnyder

Alain Person

Christian David

Arthur Du Peloux

Emeline Moubeche

Audrey Bonnelye

Pierre Dick

Moreno Eva

Arthur Du Peloux

Fault imprint in clay units: magnetic fabric, p-wave velocity, structural and mineralogical signatures

Introduction

Because of their low permeability properties coupled to a large thickness (several hundreds of meters) and a high sorption capacity, argillaceous formations have been considered as potential host rocks for nuclear waste long-term storage in several countries [START_REF] Cabrera | Caractérisation des discontinuités en milieu argileux (station expérimentale de Tournemire). Questions clés relatives à l'évaluation de la sûreté du stockage des déchets radioactifs, Rapport scientifique et technique[END_REF][START_REF] Matray | Desaturation and structure relationships around drifts excavated in the well-compacted Tournemire's argillite (Aveyron, France)[END_REF][START_REF] Sellin | The use of clay as an engineered barrier in radioactivewaste managementa review[END_REF]. However, faults and fractures in these materials are likely to impair their natural containment capabilities. For this reason, understanding the hydro-mechanical properties and the behaviour of faults in clay-rich sedimentary rocks is a topical subject for better appreciating fluid migration in deep sedimentary basins and the potential loss of integrity as geological barriers [START_REF] De Barros | Fault structure, stress, or pressure control of the seismicity in shale? Insights from a controlled experiment of fluid-induced fault reactivation[END_REF]Lefevre et al. 2017).

Anisotropy is an important characteristic that influences the clays behaviour. The properties of the argillites normally depend on the process of rock formation (deposition, compaction and diagenesis) but bedding, schistosity, and cleavage also affect the anisotropic character of the clays. Anisotropy properties of argillaceous formations can be evaluated by conducting the magnetic susceptibility p-wave velocity measurements. In clay-rich sedimentary rocks formations, the anisotropy of magnetic susceptibility (AMS) technique can be used to detect very subtle to strong tectonic deformation, even when they look undeformed at outcrop scale (Parès, 2015). On the other hand, the anisotropy of p-wave velocity (APV) is mainly controlled by the bulk elastic parameters, including the effect of pores, cracks and mineralogical content (Louis et al.2008) and is strongly influenced by deformed structures such as scaly clays, shear bands and microfolds because they are small-scale discontinuities [START_REF] Jaeggi | Tectonic structure of the "Main Fault" in the Opalinus Clay, Mont Terri rock laboratory (Switzerland)[END_REF]. Therefore, the AMS and APV are two complementary approaches for understanding the effect of failures, fractures and ductile deformation in the organization of clays. Despite the interest that these two properties have in the integral understanding of the behaviour of the rock in fracturing, there are not many studies in which these two techniques are combined. [START_REF] Louis | Anisotropy of magnetic susceptibility and P-wave velocity in core samples from the Taiwan Chelungpu-Fault Drilling Project (TCDP)[END_REF] characterized the AMS and APV in siltstones and sandstones samples in the Chelungpu fault system (Taiwan). These two techniques were also used to investigate the anisotropic behaviour of the Callovo-Oxfordian argillites [START_REF] David | Anisotropy of elastic, magnetic and microstructural properties of the Callovo-Oxfordian argillite[END_REF]) but far from a fault core zone.

In addition, the presence of faults and fractures can play an important role in controlling the migration of crustal fluids and therefore mineralogical changes are extremely sensitive to fault architecture. The intensification of the fracture network in the damage zone leads to an increased permeability with respect to the undamaged zone conducting to an authigenic growth of clays and other minerals. Conversely, the fault gouges tend to indicate a general low bulk hydraulic conductivity behaviour [START_REF] Faulkner | A review of recent developments concerning the structure, mechanics and fluid flow properties of fault zones[END_REF][START_REF] Dick | The internal architecture and permeability structures of faults in shale formations[END_REF]Lefevre et al. 2016).

Thus, it seems that a close-knit coupling exists between structure, mechanics and fluid flow properties in fault zones [START_REF] Faulkner | A review of recent developments concerning the structure, mechanics and fluid flow properties of fault zones[END_REF]. The mechanical properties of the protolith will influence the architecture of the fault zone, which will control the fluid circulation that will itself affect the mechanical properties of the fault zone.

To better understand the interrelationship between fault architecture, petrophysical properties and fluid flow circulation within a fault zone in argillaceous formation, a multidisciplinary study was made in two cores horizontally cores drilled in the Tournemire underground research laboratory (URL) in southern Aveyron (France). We used an integrated approach involving: structural analysis of plane and faults, petrophysical properties (APV and AMS) and mineralogical content.

To our knowledge this is one of the first studies to combine petrophysical and mineralogical properties of samples retrieved across a fault core zone and thus this work provides a unique opportunity to better understand the evolution of structure, petrophysical and fluid flow properties across a fault zone within clay formations. Several important questions will be addressed: Do AMS and APV properties change in phase or is their evolution different with respect to the fault? To what extent are the AMS and APV properties related to tectonic deformation? How much are these properties sensitive to localized deformation? Is it possible to observe a mineralogical change involving the presence of fluid circulation? Is there a relationship between the mineralogical changes and the petrophysical properties of the sediment?

Geological and tectonic features

The Tournemire Underground Research Laboratory (URL) is located in the South of France (Aveyron) in the western part of the so-called "Causse du Larzac" (Figure 1). The URL is operated by the French Institute for Radiological Protection and Nuclear Safety (IRSN), to study the confining properties of an argillaceous formation [START_REF] Cabrera | Structural and geophysical investigations of the EDZ (Excavation Disturbed Zone) in indurated argillaceous media: The tunnel and the galleries of the IPSN Tournemire site, France[END_REF].

The sedimentary series in this site are composed of a 250 m thick layer of clays and marls (Toarcian and Domerian) bordered by two limestone and dolomite formations (Aalenian/ Bajocian in the upper part, and Carixian/ Sinemurian in the lower part) where two aquifers are present [START_REF] Cabrera | Structural and geophysical investigations of the EDZ (Excavation Disturbed Zone) in indurated argillaceous media: The tunnel and the galleries of the IPSN Tournemire site, France[END_REF].

Fracturing evolution, associated with the regional tectonic history in the Tournemire area, was described and discussed in detail by [START_REF] Constantin | Tectonic evolution and related fracturing in the Causses Basin (Aveyron, France) : the Tournemire area example[END_REF][START_REF] Constantin | Evolution of the structural fault permeability in argillaceous rocks in a polyphased tectonic context[END_REF]. To summarize, the tectonic faults and fractures observed in the region are the results of two main tectonic phases: a first extensional tectonic phase, which comprises three episodes during the Mesozoic and a second major tectonic phase occurred during the Eocene associated to the Pyrenean compression.

The argillaceous formation of Tournemire is affected (Figure 1) by the Cernon fault (80 km long) which present a vertical and horizontal offset of several hundreds of meters and the Tournemire fault (11 km long), a local fault associated with compression tectonics [START_REF] Cabrera | Caractérisation des discontinuités en milieu argileux (station expérimentale de Tournemire). Questions clés relatives à l'évaluation de la sûreté du stockage des déchets radioactifs, Rapport scientifique et technique[END_REF].

In addition, the URL is crossed by two main parallel fault zones (F1 and F2) separated by relatively undeformed rock [START_REF] Bretaudeau | Highresolution quantitative seismic imaging of a strike-slip fault with small vertical offset in clay rocks from underground galleries: Experimental platform of Tournemire, France[END_REF] (Figure 2A). These faults of several hundred-meter length present a small vertical but a larger horizontal displacement (strike-slip faults). The main underground tunnel and several galleries intercept these faults and several boreholes have been drilled by the IRSN. Therefore, the Tournemire URL presents a unique opportunity to study in situ the evolution of the properties of clays with respect to this type of faults.

Previous studies indicate that the F1 and F2 faults show mainly a reverse left-lateral strike slip movement [START_REF] Peyaud | Mineralogical, chemical and isotopic perturbations induced in shale by fluid circulation in a fault at the Tournemire experimental site (Aveyron, France)[END_REF]. The fault's strike between N170° to N010° and 60°W to 80°W [START_REF] Lefèvre | Calcite veins as an indicator of fracture dilatancy and connectivity during strike-slip faulting in Toarcian shale (Tournemire tunnel, Southern France)[END_REF].

The sedimentary series have a sub-horizontal bedding which dips gently by 5-10° towards north but it can reach higher bedding dip values at the fault core [START_REF] Dick | The internal architecture and permeability structures of faults in shale formations[END_REF]. A dolomitic horizon marker crossed by several boreholes on either side of the fault indicated a perpendicular slip perpendicular of around 4-6 m while slip along strike is estimated to be between 15 and 30 m [START_REF] Dick | The internal architecture and permeability structures of faults in shale formations[END_REF].

Materials and strategy of sampling

We studied two boreholes (called cores in this paper): core TF1 (dip 1°, azimuth 80°N) and core ASM1 (dip 0°, azimuth 70°N). They were respectively drilled in the main tunnel in 1998 and in the "east gallery_03" in 2012 (Figure 1). ASM1 measures 101 mm in diameter and 690 cm in length and was entirely studied. TF1 measures 96 mm in diameter and 2000 cm in length but only the interval between 480 and 1432 cm and between 1578 and 1600 cm was studied. Most of the measurements were done on 2014 and 2015 and until that time, the section cores were kept vacuum sealed in foil to prevent sample alterations.

The architecture of the fault zone F2 has been already characterized in previous papers [START_REF] Dick | The internal architecture and permeability structures of faults in shale formations[END_REF], Lefevre et al. 2016). It is characterized by a central fault core and a fractured damage zone. The fault core consists in thin dark bands of centimetre thick gouge, cataclastic and brecciated rocks, as well as sub-vertical schistosity planes, folds and also lenses of less deformed rocks [START_REF] Dick | The internal architecture and permeability structures of faults in shale formations[END_REF]. The surrounding fractured damage zone is characterized by a dense network of small faults, fractures, and calcite veins. Figure 2B shows an example of the type of structures that were observed at the damage zone and the fault core zone of the TF1 fault.

The length of each zone in cores TF1 and ASM1 was determined on the basis of the description and the log carried out during the drilling operation (IRSN internal reports) (Figure 3). In both cores, the number of fractures was significantly higher in western side than in the eastern side. In addition, while the boundary between the damage and the fault core zone is sharp, the boundary between the damage and the undamaged zones was vague.

In core TF1, the western and the eastern sides of the core measure respectively 3 m and 1 m and the fault core measures 0.8m. In core ASM1, both sides have a very similar extent of around 2.5m and the fault core measures 1m (Figure 3).

After the macro fracture pattern analysis, the core sections were cut in 2x2x2cm 3 oriented cubes for the measurement of anisotropy of magnetic susceptibility (AMS, 250 samples) and the anisotropy of P-wave velocity (132 samples). Unfortunately, the section 1380-1406 was too fractured and the preparation of cubes failed.

Samples intended for bulk mineralogical content, clay mineral content and magnetic mineral analysis were obtained differently in the undamaged and the damaged zone. In the undamaged zone, a whole piece of the bulk sediment was cut and crushed in a mortar. In the damage zone and at the fault core fault, the samples were taken along the fracture planes.

About 1 mm thick layer, was scratched and crushed in a mortar into fine (bulk mineralogy, magnetic minerals) or coarse (clay minerals content) powders.

Results

Structural analysis

TF1 core is cut by numerous thin fracture planes, most of them showing a well preserved striae (Figure 4A). Most of the brittle deformation within the F2 fault zone thus corresponds to fault planes. Faulting is more intense in the damage zone (67 fractures measured in the western damage zone) than in the fault core (14 fractures measured). Due the absence of weathering, the fault surfaces are exceptionally well preserved enabling to determine the sense of slip. Seventy-six fault slip data, including fault plane and slip vector (striae), were thus collected in TF1 core, as well as five fractures with no or poorly preserved striae as shown in the rose diagram, 85% of the measured fault planes trend between N150°E and N020°E, thereby confirming the N-S orientation of the F2 fault zone and further indicating a moderate variability of the brittle fracture strike within the fault zone.

The fault population collected in TF1 core does not fit with a unique sense of shear, revealing a complex history of F2 fault. Dextral and sinistral movements, as well as normal and reverse movements exist (Figure 4A). The two last types represent 15% and 9% of the fault planes.

Strike-slip faults thus represented the majority of the measured planes (76%). Among them, dextral and sinistral faults where found in a similar proportion. The distribution direction has prominent peaks trending at N000°E -N010°E for the sinistral faults, N010°E -N020°E for the dextral and reverse faults, and N020°E -N030°E for the normal faults. This similarity in direction suggest that some faults corresponded to the reactivation of others which is further confirmed by the observation of superposed horizontal and highly dipping striae on several fault planes. As far as the fault dips are concerned, most of the N-S fault planes collected in TF1 core, whatever their movement, show value between 5° to 45° (Figure 4A). This moderate dips are rather surprising for a strike fault zone and constitute a major structural element of F2 fault.

Several mechanisms may explain these unexpected fault dips. Post-faulting local or regional tilting may have later modified the fault attitude but can be excluded here because the bedding dip does not exceed 10° and because faults have no preferred inclination polarity, dipping either to the west or to the east. An alternative scenario is that F2 fault formed as a normal fault, with eventual flattening of the brittle planes due to compaction, and was later reactivated as a strike-slip fault. Following the paleostress reconstitution of [START_REF] Constantin | Evolution of the structural fault permeability in argillaceous rocks in a polyphased tectonic context[END_REF] , the direction of extension during the Mesozoic rifting trended E-W and the direction of the Eocene compression switched from NE-SW to NW-SW direction. Such a stress succession would cause first normal slip, then dextral slip, and finally sinistral slip on F2 fault. Observations of these three types of fault movement in TF1 core suggests that F2 fault zone have suffered several periods of activities, in relation with the tectonic phases that affected the region during Meso-Cenozoic times. Alternative scenario, implying contemporaneous slips, is further developed in the discussion.

Anisotropy of magnetic susceptibility

The anisotropy of magnetic susceptibility (AMS) can be represented as a second order symmetric tensor that can be represented geometrically as an ellipsoid with three principle The range of susceptibility values found in the Tournemire sediments are comprised between 135 and 200x10 -6 SI which suggests an important contribution of paramagnetic minerals (mainly clays). Thus, the AMS ellipsoids reflect mainly the preferred orientation of clay minerals, and can be used to investigate tectonic-related fabrics [START_REF] Cifelli | The magnetic fabric in "undeformed clays": AMS and neutron texture analyses from the Rif Chain (Morocco)[END_REF][START_REF] Mattei | Magnetic fabric of weakly deformed clay-rich sediments in the Italian peninsula: Relationship with compressional and extensional tectonics[END_REF][START_REF] Parés | Evolution of magnetic fabrics during incipient deformation of mudrocks (Pyrenees, northern Spain)[END_REF]. When normalized by the density (mass magnetic susceptibility), these values ranged between 60 and 80x10 -9 m 3 /kg (Figure 5). We observed that, while in TF1 the susceptibility presented a normal distribution, a bimodal distribution was observed in ASM1. Mean values on the east (E) block of the fault were lower than in the west (W) block (Figure 5). The difference in magnetic susceptibility values between the east and west blocks could be related to the presence of calcite, a diamagnetic mineral. According to the mineralogical analyses carried out, the calcite content is higher in the eastern side. The evolution of magnetic fabric along the cores was inferred by the Jelinek diagram (T versus Pj) (Figure 5) and by the equal-area, lower-hemisphere projection of the K1 and K3 axes using the Anisoft software (Figure 6). A bedding plane dip of 10° towards the north is reported that corresponds to the measurements estimated by [START_REF] Dick | The internal architecture and permeability structures of faults in shale formations[END_REF] in the undamaged and damage zone at the west side of the fault.

With the increasing tectonic deformation, the evolution of the magnetic fabric can follow different stages [START_REF] Borradaile | Tectonic applications of magnetic susceptibility and its anisotropy[END_REF][START_REF] Parés | Evolution of magnetic fabrics during incipient deformation of mudrocks (Pyrenees, northern Spain)[END_REF][START_REF] Parés | How deformed are weakly deformed mudrocks? Insights from magnetic anisotropy[END_REF][START_REF] Parés | Sixty Years of Anisotropy of Magnetic Susceptibility in deformed sedimentary rocks[END_REF]: The type I corresponds to the sedimentary fabric and is characterized by the K3 axes perpendicular and the K1 axes are distributed on the bedding plane without a predominant magnetic lineation.

The Pj-Tdiagrams show magnetic fabrics characterized by strongly oblate shaped susceptibility ellipsoids. If the sediments are deposited under a weak current flow, the K1 axes can be aligned parallel to the paleocurrent flow azimuth.

Type II corresponds to a weak-moderate deformation with the K3 axes perpendicular to the bedding plane and the K3 axes parallel to the extensional direction. Type III, corresponding to a high deformation stage. The K3 axes show a girdle that is parallel to the maximum shortening direction and the K1 axes perpendicular to it. Furthermore, the Pj-T diagrams show a neutral ellipsoid that could be developed due to tectonic compression. The last stage of deformation (type IV) is reached when K3 is within the bedding plane and the ellipsoid is prolate.

In the undamaged zone of core TF1, the AMS ellipsoid had a predominantly oblate shape with a magnetic foliation parallel to the bedding plane (Figure 6). The shape of AMS ellipsoid at this zone had a generally oblate shape with a low degree of anisotropy (Figure 5)

and the maximum axes are relatively well clustered with a general NE-SW orientation. An oblate fabric was also observed in the damage zone. In this case the K1 axes are rotated in anti-clockwise direction. and the minimum K3 axis had a lower inclination values (75-80°)

that can be explained by the observed increase in bedding plane dip towards the fault core [START_REF] Dick | The internal architecture and permeability structures of faults in shale formations[END_REF]. Within the fault core, the AMS ellipsoids had very different signatures in sections 1353-1380 cm (west side) and 1406-1432 cm sections (east side). In section 1353-1406 cm, we can observe a type III fabric, the minimum K3 axes are distributed around a girdle from the vertical to the horizontal plane and the maximum K1 axes follows a NW-SE orientation. In addition, the Jelinek diagram (Figure 5) shows a strong decrease of the anisotropy degree (P) and a change of the magnetic shape factor (T) from high oblate (T~1)

to neutral (T~0).

The samples from section 1406-1432 section revealed also a girdle distribution of the K3 axes but the K1 axes follow a NE-SW direction. The Jelinek diagram indicates an oblate fabric.

Conversely, in ASM1, the magnetic fabric was very similar in the damage and the fault core.

It was characterized by an inclination of the K3 axes of 75-80° and an orientation of the K1 axes towards a NW-SE direction (Figure 6). We didn't observe the developpement of a type III fabric but the Jelinek diagram (Figure 5) shows a slight decrease in T from 1 to moderate values (T~0.5) indicating a weak deformation stage. In addition, we also observed higher degree of anisotropy (Pj) at the E-block. This part of the fault is also characterized by lower magnetic susceptibility (Figure 5 and Figure 7).

Anisotropy of P-wave velocity

The P-wave velocity (V P ) was measured on the cubes from both cores TF1 and ASM1 along the three axes x, y and z using ultrasonic transducers with dominant frequency 1 MHz, a

Panametrics 5058 PR pulser and a digital oscilloscope. The velocity is conventionally inferred from the travel time for a wave propagating from the transmitter to the receiver.

Previous works showed that compressional wave velocity is inversely correlated to total porosity and clay mineral contents [START_REF] Han | Effects of porosity and clay content on wave velocities in sandstones[END_REF]. In addition, in normal consolidated sediments, the Pwave velocity is faster in the bedding plane than perpendicular to it.

The anisotropy of P-wave velocity (APV) is defined as the difference between maximum and minimum velocities normalized to the average velocity:

APV(%)=100*(Vp max-Vp min)/Vp mean, where Vpmin is the P-wave velocity in the Z direction (VPz), perpendicular to bedding, and Vpmean the average velocity in the bedding plane (Vpx+Vpy)/2. The Vp mean, Vpz and the APV (%) are shown in Figure 7. The P-wave velocities in the bedding plane were much higher than the vertical one and lead to large APV values of around 50% in the undeformed zone to around 40% in the damage zone. In core TF1, we observed a drastic change within the fault core with a strong dispersion of APV values down to negative values, reflecting elastic anisotropy reversal. This is mainly due to the increase in P-wave velocity along the vertical axis and a decrease in the horizontal axes.

Such a change is not observed in the fault core of ASM1. We can observe that in general the changes in the P-wave velocity and the APV are in phase with changes the AMS properties (Figure 7).

Clay minerals assemblages and bulk mineralogy

Bulk mineralogy and clay mineral assemblages were identified using X-ray diffraction (XRD). It should be remembered that the samples for the XRD and MEB measurements come from a block of total sediment in the undamaged zone and from the first millimetres along the slip planes from fractures in the damage zone and in the fault core.

XRD samples were prepared according to the powder diffraction method [START_REF] Brindley | Crystal structures of clay minerals and their X-ray identification[END_REF]. For clay mineral preparation, we followed the analytical procedure of [START_REF] Holtzapffel | Les minéraux argileux : préparation, analyse diffractométrique et détermination[END_REF]. X-ray diffractograms were obtained using a D2 Brucker diffractometer equipped with a Lynx Eye detector, with CuKα radiation and NI filter. Measurement parameters were as follows: 2.5 to 35° (2θ), in steps of 0.02° each 0.2 s. The identification of clay minerals was made according to the position of the (001) reflections of the X-ray diffractograms [START_REF] Brindley | Crystal structures of clay minerals and their X-ray identification[END_REF][START_REF] Moore | X-ray diffraction and the identification and analysis of clay minerals[END_REF]) using a MacDiff software.

XRD analysis was complemented by the analysis of some bulk samples from the undamaged host rock and fractures planes within the damage zone in a scanning electron microscopy coupled to an energy-dispersive X-ray spectrometry system (SEM-EDS). We used a ZEISS supra V5 instrument provided with a BRUKER detector in energy dispersive spectroscopy (EDS) at the Petrology, Geochemistry, Magmatic Mineralogy Laboratory (PG2M) ISTEP of the University Pierre and Marie Curie (Paris).

The major mineral assemblages identified by XRD analysis were clays (18-48%, mean 37%), quartz (21-48%, mean 37%) and calcite (9-61%; mean 27%). Pyrite and siderite were found as secondary minerals with concentrations below 5% and they were not quantified. We observed small quantities of ankerite and ferroan dolomite both in the undeformed host rock sediment and in the damage zone by SEM-EDS (Figure 9). This finding is in agreement with previous studies [START_REF] Lerouge | Diagenetic carbonates in clay-rich marine formations[END_REF][START_REF] Peyaud | Mineralogical, chemical and isotopic perturbations induced in shale by fluid circulation in a fault at the Tournemire experimental site (Aveyron, France)[END_REF]. We must note that core TF1 was drilled in 1999 and that some mineralogical alteration could have occurred. However, its composition is similar to the one in core ASM1, collected 14 years later indicating that postdrilling alteration didn't change the original mineralogical composition.

The calcite content in the western undamaged zone is very low (~15%) (Figure 8) and corresponds to the average calcite content estimated in the Tournemire site [START_REF] Boisson | Caractérisation des propriétés de transferts au laboratoire des argilites et marnes toarciennes de Tournemire[END_REF]). As we approach the centre of the fault, the calcite content from slip planes showed very scattered values (10-60% in core TF1 and 9-51% in core ASM1) indicating that some fractures are filled with calcite cement and some other are not. This calcification took place during two episodes of compressive deformation during the Pyrenean orogeny [START_REF] Lefèvre | Calcite veins as an indicator of fracture dilatancy and connectivity during strike-slip faulting in Toarcian shale (Tournemire tunnel, Southern France)[END_REF].

Additionally, we also observed higher CaCO 3 content in the eastern compartment (30-35%) (end of damage zone and undeformed zone), reaching 30-35% (Figure 8). We interpret that this increase is not caused by precipitations from tectonically induced fluids (macroscopically, we found no veins in this core section/ these samples) but that it rather indicates a different background signature for calcite content eastwards of the fault zone.

Average proportions of the various clay species in TF1 and ASM1 cores are very similar (Figure 8). The clay mineral assemblage in core TF1consisted of illite (35 to 45%, near 40% on average), kaolinite (25 to 40%, 33% on average), illite-smectite mixed layers (I/S, 15 to 25%, 19% in average) and chlorite (5 to 10%, 9% in average).

The clay mineral assemblage in core ASM1 consisted of illite (5 to 40%, 36% on average), kaolinite (30 to around 40%, 34% on average), I-S (15 to 35%, 21% on average) and chlorite

(5 to 10%, 9% on average).

Clay mineral proportions in TF1 and ASM1 cores were also very similar to the ones previously determined at Tournemire site [START_REF] Savoye | Origin of pore-water isotopic anomalies near fractures in the tournemire shales[END_REF][START_REF] Dick | The internal architecture and permeability structures of faults in shale formations[END_REF]). However, while [START_REF] Savoye | Origin of pore-water isotopic anomalies near fractures in the tournemire shales[END_REF] didn't show any significant change in the clay assemblage as function of the distance to the fault plane, our results revealed a weak but noticeable increase in kaolinite (Figure 8). The maximal difference in kaolinite content between the undamaged and the fault core zones is 14% in core ASM1 and 9%, in core TF1, slightly higher than the analytical error (±5%). We interpret this mineralogical evolution as significant as the same mineralogical trend is observed in both ASM1 and TF1 cores toward the deformed zones. A previous work indicates that the proportion ok kaolinite, can be even higher (63-72%) when we look at the fault gouge [START_REF] Dick | The internal architecture and permeability structures of faults in shale formations[END_REF]. Such a kaolinite increase toward the fault core has been considered to reflect a high degree of fluid-rock interaction (Rossetti et al. 2010).

An additional mineralogical change recorded in our data set ASM1 (Figure 8) is characterized by a large increase in I/S proportions correlated with a decrease in illite from the fault core boundary (I/S 16% and illite 42%) to the eastern undamaged zone (I/S 35% and illite 23%). This sharp mineralogical change may suggest a vertical offset that can be more important in ASM1. After the difference in the dolomitic horizon marker observed on either side of the fault, suggest that the amount of the perpendicular slip is of around 4-6 m (Dick et al. 2016).

Magnetic mineralogy based on low temperature SIRM

Low-temperature isothermal remanence magnetizations (LT-SIRM) were made in core TF1 using a Quantum Design Magnetic Property Measurement System (MPMS2) SQUID magnetometer at the IPGP Paleomagnetic laboratory in Paris (France). The powders used for these analysis were collected on the fault planes and fractures at different points of core TF1

with the exception of the sample in the undeformed zone which comes from the bulk sediment.

Two different LT-SIRM experiments were performed. First, samples were cooled in zero field to 10 K, and a magnetic field of 2.5T was given at 10 K. The SIRM was measured in approximately zero field at 5-K intervals up to 300 K (zero-field cool curves, ZFC). In the second experiment, the samples were cooled in a high magnetic field of 2.5 T and then the SIRM was measured again every 5K up to 300 K (field-cool curves, FC).

We observed two transitions in the ZFC and FC curves: A small drop at 120 K, in the majority of the samples (sometimes only observed in the first derivate curves) which corresponds to temperature of the magnetite Verwey transition, Tv [START_REF] Verwey | Electronic Conduction of Magnetite (Fe3O4) and its Transition Point at Low Temperatures[END_REF] (Figure 10). The elevation of FC SIRM over ZFC SIRM and the small loss in remanence across the Tv warming indicates a very fine grain-size, lower than 0.1µm which correspond to singledomain (SD) magnetite [START_REF] Özdemir | Changes in remanence, coercivity and domain state at low temperature in magnetite[END_REF][START_REF] Smirnov | Grain size dependence of low-temperature remanent magnetization in natural and synthetic magnetite: Experimental study[END_REF]. In addition, a linear decrease trend observed in the FC and ZFC SIRM curves above the Verwey transition (T > 120 K)

indicates the presence of superparamagnetic (SP) grains. The frontier between the SP and SD grains depends on the shape of magnetite but for detrital magnetite with isotropic grains, the limit is at around 0.017 µm [START_REF] Muxworthy | Critical superparamagnetic/single-domain grain sizes in interacting magnetite particles: implications for magnetosome crystals[END_REF].

In addition, magnetite is characterized to have a high remanence at room temperature and we used the SIRM at 300K (SIRM 300K ) to evaluate the evolution of magnetite content along the TF1.

A second large drop in remanence was observed around 35K and it was related to the Neel transition of siderite. Siderite is also characterized by a FC remanence below 40K much larger than the ZFC remanence [START_REF] Housen | Low-temperature magnetic properties of siderite and magnetite in marine sediments[END_REF].

Two other minerals have magnetic transitions around 35K: rhodochrosite (MnCO 3 ) whose Neel temperature is 32K [START_REF] Housen | Low-temperature magnetic properties of siderite and magnetite in marine sediments[END_REF] and pyrrhotite (Fe 1-x S) which has a magnetic transition at 35 K [START_REF] Dekkers | Grain-size dependence of the magnetic behavior of pyrrhotite during its low-temperature transition at 34 K[END_REF][START_REF] Rochette | Magnetic transition at 30-34 K in pyrrhotite-insigth into a widespread occurence of this mineral in rocks[END_REF] The magnetic mineralogy evolution from the undeformed host rock to the fault core was studied using the SIRM 300K used as an estimate of fine-grained magnetite and FC/ZFC ratio at 10 K, used to indicate the occurrence of siderite. For pure siderite, the difference is higher than 10 [START_REF] Frederichs | Towards the identification of siderite, rhodochrosite, and vivianite in sediments by their low-temperature magnetic properties[END_REF][START_REF] Housen | Low-temperature magnetic properties of siderite and magnetite in marine sediments[END_REF].

We observed a sharp increase in the SIRM 300K values in the damage zone (Figure 10) in phase with a decrease in the FC/ZFC 10K ratio. These changes were interpreted as a partial oxidation of siderite to magnetite due to the presence of fluid circulations [START_REF] Ellwood | Anomalous magnetic properties in rocks containing the mineral siderite: Paleomagnetic implications[END_REF].

Discussion

The multi-proxies approach used in this study confirmed the commonly accepted zooning along a fault zone: an undamaged area characterized by samples with little or no fracturing, a damage zone with a large number of calcified fractures and a fault core characterized by high deformation and no apparent bedding. Based on this finding, it is tempting to assume a fault zone evolution where a tectonic rework of the host rock lead to the formation of gouge, which we assume to accommodate the largest part of the faults' offset [START_REF] Dick | The internal architecture and permeability structures of faults in shale formations[END_REF][START_REF] Laurich | Deformation mechanisms and evolution of the microstructure of gouge in the Main Fault in Opalinus Clay in the Mont Terri rock laboratory (CH)[END_REF].

Our work indicates three main findings:

-The tectonic analysis of fractures indicated that the F2 fault underwent a polyphased evolution with reactivations under different tectonic regimes. Two strike-slip regimes, with approximate NE-SW and NW-SE compression, and a normal regime with a WNW-ESE direction of extension are inferred from the orientation of the fault-slips. Most of the N-S fault planes collected in TF1 core, whatever their movement, show dip values between 5° to 45°. This finding constitutes a major structural characteristic of F2 fault.

-Second, our study demonstrated that the anisotropy of P-wave velocity (APV) and of susceptibility (AMS) follows the fault zoning in agreement with the results obtained by [START_REF] Bonnelye | Etude des propriétés physiques et mécaniques des argilites : de la déformation en laboratoire aux failles naturelles[END_REF]. We observed the loss of the primary sedimentary fabric in the damage zone and an important decrease in the anisotropy of magnetic susceptibility (AMS) and the P-wave velocity (APV) within the fault core, especially in TF1 core.

-Third, we identified evidence of neo-formation of fine-grained magnetite (only measured in core TF1), precipitation of millimetric to infra-millimetric calcite veins and a potential neoformation of kaolinite within the fault planes in the damage and the core fault.

Polyphased fault zone development

Fault data collected along the F2 fault evidenced the complex nature of the damage zones of this fault. The damage zones include individual fractures with horizontal and dip slip movements, sometimes in opposite directions. The NS+-20 sinistral slip surfaces observed in this study are in agreement with a left lateral movement observed on other segments of F2 fault in the Tournemire URL [START_REF] Constantin | Tectonic evolution and related fracturing in the Causses Basin (Aveyron, France) : the Tournemire area example[END_REF][START_REF] Lefèvre | Calcite veins as an indicator of fracture dilatancy and connectivity during strike-slip faulting in Toarcian shale (Tournemire tunnel, Southern France)[END_REF].

Sinistral and reverse faults in TF1 both indicate a local N130°E  H (maximal horizontal stress) and a N040°E h (minimal horizontal stress), in agreement with the N130°E compression inferred from fault-slip data inversion close to the northern continuation of F2 in the western gallery [START_REF] Constantin | Tectonic evolution and related fracturing in the Causses Basin (Aveyron, France) : the Tournemire area example[END_REF]. Dextral faults result from a N45°E  H and N135°E  h . A rather similar WNW-ESE  h is inferred from the normal faults and is associated with a normal stress regime.

We propose two different scenarios that may explain this tectonic complexity (Figure 11). We consider (1) a long lived scenario, with fault activity starting in Mesozoic time and regional stresses changing through time (Figure 11A), and (2) a short scenario, in which the fault internal architecture evolved progressively during the displacement accumulation (Figure 11B). Scenarios between these two end member models may further exist. In scenario 1 (Figure 11A), the F2 fault initiated as a normal fault during Mesozoic time under a WNW-ESE extension. As such an extension regime is active during the Toarcian [START_REF] Constantin | Tectonic evolution and related fracturing in the Causses Basin (Aveyron, France) : the Tournemire area example[END_REF], normal faults in Tournemire may predated sedimentary compaction. Subsequent, the F2 fault registered a first dextral slip and a later sinistral one in Eocene time during the Pyrenean orogeny. These two opposite slips result from the counter clockwise stress rotation of the Pyrenean compression, from a N20°E to a N160°E direction [START_REF] Constantin | Tectonic evolution and related fracturing in the Causses Basin (Aveyron, France) : the Tournemire area example[END_REF].

In scenario 2 (Figure 11B), the F2 fault observed in the gallery section is a linking fault in an extensive relay between two NNW-SSE dextral faults formed during the Pyrenean orogeny.

The normal faults recognized in the TF1 core formed when these two dextral segments were linked (Figure 11B1). Such a development is confirmed in analogue models [START_REF] Rahe | Physical analog modeling of pull-apart basin evolution[END_REF]. In the latest stage, the strike-slip fault cut the relay (Figure 11B2) and formed the F2 fault. The relay zone was later reactivated as a sinistral strike-slip fault in response to the counter clockwise rotation of the Pyrenean compression (Figure 11B3).

Both scenarios are in agreement with the sinistral reactivation of previously dextral faults observed by Lefevre et al. (2016) along another section of F2 fault. Interestingly, the earlier dextral slips are more preserved in the studied portion of the fault zone (50% of the measured fault surfaces in TF1 core). The left lateral reactivation of the fault requires that N160°E compression was deflected close to the fault in order to have a sufficient shear stress.

Accordingly, the local stress tensor obtained by [START_REF] Constantin | Tectonic evolution and related fracturing in the Causses Basin (Aveyron, France) : the Tournemire area example[END_REF] in the Tournemire gallery indicates that the compression adopted here a N130°E strike.

Fault zoning and longitudinal heterogeneity of anisotropy of p-wave velocity and magnetics susceptibility

We found an excellent correlation between the anisotropy of magnetic susceptibility (AMS)

and the P-wave velocity (APV) in cores TF1 and ASM1. In the case of argillites, the susceptibility and P-wave velocity are mainly controlled by the preferred orientation of the clay sheets and explains while the magnetic and elastic results converge [START_REF] David | Anisotropy of elastic, magnetic and microstructural properties of the Callovo-Oxfordian argillite[END_REF].

The results of AMS in ASM1 and TF1 cores indicate that there is a dominance of type II and type III fabrics with oblate to neutral ellipsoids. A magnetic lineation trending NW-SE was observed in the damage zone and in majority of the samples of the fault core zones. Following [START_REF] Mattei | Magnetic fabric of weakly deformed clay-rich sediments in the Italian peninsula: Relationship with compressional and extensional tectonics[END_REF] and [START_REF] Parés | How deformed are weakly deformed mudrocks? Insights from magnetic anisotropy[END_REF], this magnetic lineation is usually developed parallel to the extension direction. The observed K1 preferred orientation can result from strike-slip regime with  1 trending NE-SW. Such stress state is in agreement with the numerous dextral slip measured in the fault zone. The AMS fabric was thus acquired during the dextral movement along F2 fault during the Pyrenean orogeny. The normal movement, whatever their age, may also have produced the magnetic fabric, but its contribution is thought to have been less important.

Interestingly, the magnetic fabric in section 1402-1432 in the TF1 core showed a different pattern. In this case, the K3 axes formed a girdle in the NW-SE direction and the K1 axes were preferably oriented in the NE-SW direction. Several mechanisms could explain this change. It could indicate that the eastern side of the fault imprints the deformation due to sinistral reactivation of previously dextral faults. This fabric was also observed in this section. It would be necessary to study other parts of the fault to affirm or decline this hypothesis. On the other hand, AMS tectonic fabric is very sensitive to localized deformation zones such as gouge or deformation bands that are very abundant in this part of the fault. The impact of these deformations on the AMS fabric is an issue that needs to be further explored [START_REF] Parés | Sixty Years of Anisotropy of Magnetic Susceptibility in deformed sedimentary rocks[END_REF].

Similarly, a complete reversal of elastic anisotropy (APV <0) was also observed within the fault core. Recent experiments on Tournemire shale [START_REF] Bonnelye | Elastic wave velocity evolution of shales deformed under uppermost-crustal conditions[END_REF] and creep tests [START_REF] Geng | Elastic anisotropy reversal during brittle creep in shale[END_REF] show an inversion of elastic anisotropy in samples deformed perpendicularly to bed changes due to mineral rotation near fault zones. A tectonic fabric is also detected in a borehole from the Mont Terri rock laboratory thanks to a decrease in the pwave velocity [START_REF] Jaeggi | Tectonic structure of the "Main Fault" in the Opalinus Clay, Mont Terri rock laboratory (Switzerland)[END_REF].

Thus, the petrophysical results observed within the fault core of TF1 were interpreted as the result of an intense damage process. Here, the bed is only preserved in centimetre -thick zones (Lefevre et al. 2016) which may explain why the % APV and K3 axes showed such scattered values (Figure 7), these parameters being particularly sensitive to bedding orientation.

Another important result is that the preferred orientation of the K1 axes in the undamaged zone (only available in TF1) was in the NE-SW direction. One possible mechanism is that the alignment of the K1 axes in the undamaged zone is due to a deposition under the influence of a bottom current. In this case, the direction of the K3 axes is perpendicular to the bedding plane and the K1 axes are parallel to the current direction [START_REF] Beckers | Influence of bottom currents on the sedimentary processes at the western tip of the Gulf of Corinth, Greece[END_REF][START_REF] Singsoupho | Palaeocurrent direction estimated in Mesozoic redbeds of the Khorat Plateau, Lao PDR, Indochina Block using anisotropy of magnetic susceptibility[END_REF]Parès et al. 2007). Further studies will be necessary to confirm this hypothesis.

In ASM1, the AMS and elastic anisotropy didn't show such drastic changes in the fault core.

The AMS analysis indicated a type II tectonic fabric [START_REF] Parés | Sixty Years of Anisotropy of Magnetic Susceptibility in deformed sedimentary rocks[END_REF] and a slight decrease in the p-wave velocity at 300 cm within the fault core. In the fault core, the sediment is characterized by small deformed volumes, shear zones, shear bands and gouge intercalated with apparently undeformed sediment [START_REF] Jaeggi | Tectonic structure of the "Main Fault" in the Opalinus Clay, Mont Terri rock laboratory (Switzerland)[END_REF][START_REF] Dick | The internal architecture and permeability structures of faults in shale formations[END_REF]Lefevre et al. 2016). Despite the care we took in the sampling procedure, the number of sampled cubes were limited by the large number of fractures and it cannot be totally ruled out that the deformation of these relatively small volumes zones were less reflected in ASM1 because they couldn't be sampled. The magnetic fabric in the damage and the fault core in ASM1

show however a similar a preferred orientation of the K1 awes in the NW-SE direction in agreement with the dextral strike-slip regime.

Finally, one of the major characteristic found in core ASM1 was the dissymmetry observed in P-wave velocity, magnetic susceptibility and Pj between the western and eastern parts of the fault. The observed differences seem to be related to a mineralogical change and not to strain and they were not observed in TF1. Mineralogical analysis indicates higher calcite content in the eastern side and an increase in illite-smectite mixed layers towards the east. As the core was drilled horizontal to the bedding plane, the stratigraphic level (and therefore the mineralogical composition), should be the same on both sides of the fault, unless there is a vertical offset. This vertical offset was observed by [START_REF] Dick | The internal architecture and permeability structures of faults in shale formations[END_REF]. We cannot however quantify it.

Fluid circulation in the fault zone

The XRD analyses showed evidence of precipitation of millimetric to infra-millimetric calcite veins and a potential neoformation of kaolinite along the fault planes in the damage zone/fault core zones. This was already observed in a previous study of Tournemire shales which indicated that maximum kaolinite content was observed in the fault gouge [START_REF] Dick | The internal architecture and permeability structures of faults in shale formations[END_REF].

The increase of kaolinite abundance towards the fault core suggests a high degree of fluidrock interaction [START_REF] Rossetti | Fluid flow within the damage zone of the Boccheggiano extensional fault (Larderello-Travale geothermal field, central Italy): Structures, alteration and implications for hydrothermal mineralization in extensional settings[END_REF].

In addition, the LT-SIRM measurements enabled the identification of small amount of siderite in our samples. Siderite was also identified by XRD but it was not quantified by XRD because its concentration was less than 5%. However, magnetic measurements like SIRM are able to detect ferrimagnetic minerals that are undetected by traditional mineralogical analyses such as XRD, even when the concentrations are as low as 1 ppm [START_REF] Lagroix | A New Tool for Separating the Magnetic Mineralogy of Complex Mineral Assemblages from Low Temperature Magnetic Behavior[END_REF].

Besides, the SIRM 300k and FC/ZRC record indicated an increase in fine-grained magnetite and a decrease in siderite towards the damage/fault core boundary (Figure 9).

In Tournemire shales, the fluid circulation mainly occurs between the fault core and the damage zone [START_REF] Guglielmi | In situ observations on the coupling between hydraulic diffusivity and displacements during fault reactivation in shales[END_REF][START_REF] Lefèvre | Calcite veins as an indicator of fracture dilatancy and connectivity during strike-slip faulting in Toarcian shale (Tournemire tunnel, Southern France)[END_REF] in which the porosity reaches maximum values [START_REF] Dick | The internal architecture and permeability structures of faults in shale formations[END_REF]. Therefore, the transformation of siderite into magnetite seems to be a good index of paleo-fluid circulation and indicates a change in redox conditions related with episodic connections with less reduced fluids from the overlying aquifer [START_REF] Peyaud | Mineralogical, chemical and isotopic perturbations induced in shale by fluid circulation in a fault at the Tournemire experimental site (Aveyron, France)[END_REF]. We suggest however the prevalence of reducing conditions, as at higher oxidation rates, the final product of siderite will be hematite instead of magnetite [START_REF] Ellwood | Anomalous magnetic properties in rocks containing the mineral siderite: Paleomagnetic implications[END_REF]. The presence of minerals like ankerite, calcite, ferroan dolomite and framboidal pyrite observed in the SEM (Figure 8) supports also the prevalence of reducing conditions and were also observed by [START_REF] Peyaud | Mineralogical, chemical and isotopic perturbations induced in shale by fluid circulation in a fault at the Tournemire experimental site (Aveyron, France)[END_REF] and [START_REF] Lefèvre | Calcite veins as an indicator of fracture dilatancy and connectivity during strike-slip faulting in Toarcian shale (Tournemire tunnel, Southern France)[END_REF].

Conclusions

The integrated study of the N-S F2 fault cutting the Toarcian shale formation from the Tournemire Underground Research Laboratory (URL) made it possible to highlight the existing relationships between the fault architecture and its tectonic history with the variations of petrophysical properties (anisotropy of magnetic susceptibility and P-wave velocity) and mineralogy (CaCO 3 , clays and magnetic mineralogy).

The fault surfaces within the fault zone pointed out that the F2 fault suffered a polyphased tectonic history, with successive reactivations under normal and strike-slip regimes. Notably, slips on fault surface within the damage zone occurred mostly on shallow or very shallow (5° to 45°) dipping N-S fault planes.

Our study demonstrated that the anisotropy of P-wave velocity (APV) and of susceptibility (AMS) are well correlated and they follow the fault zoning. In the damage zone, changes in the ASM and APV properties were mainly related to an increase in the bedding plane dip.

Within the fault core zone, a strong anisotropy reversal was observed which indicated a strong sediment deformation and suggest clay platelets rotation near failure. These deformations are however located in very small volumes (deformed bands and gouge). [START_REF] Mennessier | Carte géologique de la France[END_REF]. (B) Geological cross-section of the Tournemire experimental station [START_REF] Cabrera | Caractérisation des discontinuités en milieu argileux (station expérimentale de Tournemire). Questions clés relatives à l'évaluation de la sûreté du stockage des déchets radioactifs, Rapport scientifique et technique[END_REF].

  susceptibility axes: maximum (K1), intermediate (K2) and minimum (K3) susceptibility with mean magnetic susceptibility Km=(K1+K2+K3)/3. The AMS was estimated using a Geofyzika KLY-3 KappaBridge at the Paleomagnetism laboratory in the Institut de Physique du Globe de Paris (IPGP) (France).

Jelinek ( 1981 )

 1981 introduced different parameters to interpret the degree and shape of the anisotropy of magnetic susceptibility: lineation (L=K1/K2); foliation (F=K2/K3). The shape parameter (T) combines the lineation and foliation parameters (T = ((ln L -ln F)/(ln L + ln F)). Positive T-values imply an oblate ellipsoid, negative values imply prolate ellipsoid, and T-values close to zero imply neutral ellipsoids. Finally, Pj is the corrected degree of anisotropy (Pj= exp( 2[(ln(K1/km) 2 +ln(K2/km) 2 +ln(K3/km) 2 ]). Higher Pj values imply increasing strength of the ellipsoid shape.

FranceFigure 1 :

 1 Figure 1: (A) Location of the Tournemire site: simplified geological map after[START_REF] Mennessier | Carte géologique de la France[END_REF]. (B) Geological cross-section of the Tournemire experimental station[START_REF] Cabrera | Caractérisation des discontinuités en milieu argileux (station expérimentale de Tournemire). Questions clés relatives à l'évaluation de la sûreté du stockage des déchets radioactifs, Rapport scientifique et technique[END_REF].

Figure

  Figure (with caption below and on the same page)

FigureFigure 3 :

 3 Figure (with caption below and on the same page)

FigureFigure 4 :

 4 Figure 3: Logs of ASM1 (A) (société Martinez and Pierre Dick, personal communication) and TF1 (B)(Cofor -GeoSonic France for the IRSN) made during the drilling process. The spatial distribution of the undamaged, damage and fault core zones was determined on observation criteria during the core logging. The number cubes for ASM and (APV) are indicated. The black and white circles correspond respectively to the near shear points in where the powder samples were collected for the XRD and low-T SIRM measurements. Figure (with caption below and on the same page)

FigureFigure 5 :

 5 Figure (with caption below and on the same page)

Figure

  Figure (with caption below and on the same page)

Figure 6 :

 6 Figure6: Equal-area lower-hemisphere stereograms in geographical coordinate system of the principal AMS axes for the undamaged zone, damage zone and fault core zone of cores TF1 and ASM1. White circles represent K3 axes, triangles and black squares represent K1 axes. The black line represents a bedding plane oriented to the north with a 10°dip which corresponding to the values measured in the undamaged zone.

Figure 7 :

 7 Figure (with caption below and on the same page)

FigureFigure 8 :

 8 Figure (with caption below and on the same page)

Figure (with captionFigure 9

 9 Figure (with caption below and on the same page)

Figure

  Figure 9: (A) SEM image of a broken surface in core TF1, section 980-1000 (undamaged zone). (B) Si, S, Fe, Ca elements mapping of the same sample. Silicates, calcites, pyrites and small quantities of iron carbonates were observed. (C) EDS spectrum of iron reach carbonate observed in SEM image.Figure (with caption below and on the same page)

Figure 10 :

 10 Figure 10: LT-SIRM warming curves measured in powders from the fracture planes in the damage zone (sample 1180) and fault core zone (sample 1353) in core TF1. FC (Field cool) and ZFC (Zero field curve) correspond respectively to samples cooled at zero field and a field of 2.5T. (B) Evolution of the SIRM at room temperature (SIRM300K) and FC/ZFC ratio at 10K across the undamaged, damage and the fault core zones.

FigureFigure 11 :

 11 Figure (with caption below and on the same page)

  Figure (with caption below and on the same page)

  now known as Benus transition[START_REF] Rochette | The Low-Temperature Magnetic Transition of Monoclinic Pyrrhotite[END_REF].[START_REF] Aubourg | Toward a new < 250 °C pyrrhotite-magnetite geothermometer for claystones[END_REF] as well as[START_REF] Aubourg | Imprinting chemical remanent magnetization in claystones at 95 °C[END_REF] attributed the drop in the ZFC/FC curves in samples from Mont Terri Lower Dogger claystones to a

combination of paramagnetic input (clay minerals), pyrrhotite transition and the effect of superparamagnetic grains (SP). However, in Tournemire site, pyrrhotite wasn't found neither by XRD analysis nor SEM/EDS observations while siderite can be detected by XRD and have concentrations lower than 5%.

Acknowledgment

This work was funded by a collaboration contract between the IRSN, the CNRS and the MNHN (number DAFCJ/SL/BL-11831). We are grateful to France Lagroix for allowing us the access to the Paleomagnetism laboratory of the IPGP for the susceptibility measurements and for her support given during the low-T SIRM measurements. We are grateful and to Omar Boudouma for his support given during the SEM sessions and for Frederic Delbès in helping for clay mineral analyses. We thank Ben Laurich and the anonymous reviewer for their constructive criticism of the paper.

The observed NW-SE magnetic lineation in both cores ASM1 and TF1 was acquired during the first dextral slip during the Pyrenean orogeny.

The magnetite and potential kaolinite neoformation attest of fluid circulation at the vicinity of the damage/fault core boundary and suggest that the transformation of siderite in magnetite may be a good index of paleo-fluid circulation.

The variability of the response observed in the two studied cores suggests the need to study the characteristics of the fault both across its different architectural zones and along the main fracture plan.