Seismic anisotropy of Opalinus Clay: tomographic investigations using the infrastructure of an underground rock laboratory (URL)

Seismic anisotropy and attenuation make claystone formations difficult to characterize. On the other hand, in many geotechnical environments, precise knowledge of structure and elastic properties of clay formations is needed. In crystalline and rock salt underground structures, high-resolution seismic tomography and reflection imaging have proven a useful tool for structural and mechanical characterization at the scale of underground infrastructure (several deca- to hundreds of meters). This study investigates the applicability of seismic tomography for the characterization of claystone formations from an underground rock laboratory under challenging on-site conditions including anisotropy, strong attenuation and restricted acquisition geometry. The seismic tomographic survey was part of a pilot experiment in the Opalinus Clay of the Mont Terri Rock Laboratory, using 3-component geophones and rock anchors, which are installed 2 m within the rock on two levels, thus suppressing effects caused by the excavation damage zone. As a source, a pneumatic impact source was used. The survey covers two different facies types (shaly and carbonate-rich sandy), for which the elliptical anisotropy is calculated for direct ray paths by fitting an ellipse to the separated data for each facies. The tomographic inversion was done with a code providing a good grid control and enabling to take the seismic anisotropy into account. A-priori anisotropy can be attributed to the grid points, taking various facies types or other heterogeneities into account. Tomographic results, compared to computations using an isotropic velocity model, show that results are significantly enhanced by considering the anisotropy and demonstrate the ability of the approach to characterize heterogeneities of geological structures between the galleries of the rock laboratory.


Introduction
In the development and usage of underground infrastructure, claystones are particularly challenging, but they also have relevant properties which make them useful for a wide range of geoengineering applications. Their extremely low permeability has enabled the generation of hydrocarbon reservoirs in structural traps (e.g., Grunau 1987), and the same property is used in the context of geological CO 2 storage in deep saline aquifers, located underneath one or (usually) more layers of argillaceous rock formations (e.g., Espinoza and Santamarina 2017). Their retention capacity, and their ability to swell are considered favourable properties in nuclear waste disposal (Mallants et al. 2001;Delage et al. 2010), whereas these properties, especially swelling, may pose severe challenges in the construction of underground infrastructure (e.g., Einstein 1996).
Therefore, a detailed characterization of clay formations is of high importance in many applications. Mechanical and hydraulic properties, especially potential heterogeneities, have significant impact on the assessment of short-and long-term behaviour of these formations under varying environmental conditions.
The clay mineral composition has in some instances a large impact on the physical properties with just a small amount of certain clay minerals (Murray 2006). Furthermore a complex anisotropy (structural anisotropy) results from interferences of different components of the rock structure (e.g. bedding planes, clay mineral-rich layers, sandy lenses, and cemented matrix). An example are laminated materials like Opalinus Clay (Siegesmund et al. 2014).
Underground geophysical exploration using seismic sources and receivers has been applied successfully in various geological settings including crystalline or rock salt formations (e.g. Richter et al. 2018), but little experience is available for underground rock laboratory (URL) scale characterization in claystone, mostly due to the above mentioned particular challenges posed by this rock type.
The rock laboratory of the Mont Terri Project (www. mont-terri. ch, Bossart et al. 2017), located near St-Ursanne in the canton of Jura (Switzerland), is offering excellent experimental conditions for the acquisition of geophysical data in for characterization at a large range of scales. A wide range of experiments have been performed since the foundation to research on the hydrogeological, geochemical, and geotechnical properties to characterize the Opalinus Clay formation. Bossart et al. (2017) give a thematic overview for the large variety of experiments held from 1996 to 2016. Among those are many seismic investigations, e.g., to characterize changes in rock properties during excavation (e.g. Le Gonidec et al. 2012;Schuster et al. 2019), to study resulting changes in form of the excavation damage zone (EDZ, e.g. Schuster et al. 2001;Nicollin et al. 2008;Schuster et al. 2017) and experiments for seismic monitoring (e.g. Manukyan et al. 2012;Zappone et al. 2021). Another research aspect is the general characterization of seismic velocities and anisotropy of the Opalinus Clay and its facies in situ and in the laboratory (e.g. Popp and Salzer 2007;Sarout et al. 2014;Schuster et al. 2017). Most of these experiments put a stronger focus onto a very high resolution but small scale characterization. When larger scale structures were illuminated, usually rather sparse acquisition geometries were applied providing average effective properties over a large distance.
The aim of this study was to perform a tomographic inversion of this relatively sparse seismic data set and to investigate the feasibility of seismic characterization of the different Opalinus Clay facies types. Heterogeneities within the shaly and the carbonate-rich sandy facies were not expected to be resolved in detail, but the investigation should demonstrate that the two facies types could be discriminated from each other by tomographic traveltime inversion and the undisturbed rock masses could be characterized in terms of compressional and shear wave propagation velocities.

Geological setting
The Mont Terri anticline is located in the Folded Jura, which is part of the Jura mountain belt, an arcuate foldand-thrust belt Jaeggi et al. 2017). The URL is located in the southern limb of the anticline embedded in the Opalinus Clay formation, an overconsolidated claystone, with an overburden of 250-320 m . The bedding of the formation dips in SSE direction and the dip decreases from around 50 • to 30 • from the northern to the southern part of the laboratory (Nussbaum et al. 2011).
The EDZ is a well characterized feature along the tunnel walls. With a depth of 1-2 m in gallery 08 (Lanyon et al. 2014) we assume a similar depth for niche 8 which has the same orientation as gallery 08 perpendicular to the bedding. For the niche TT to gallery 18 we have a orientation parallel to the bedding but for the carbonaterich sandy facies no macroscopic EDZ fractures are identified (Lanyon et al. 2014).

Seismic survey
In January 2019, we carried out a seismic pilot survey at the Mont Terri Laboratory that covered the galleries 04 and 08 as well as the niches TT, 7, and 8 ( Fig. 1B, red rectangle). The pilot survey aimed at the evaluation of the applicability and performance of different sources in clay formations (Wawerzinek et al. 2021). The second major aspect of the survey and focus of this study was the characterization of claystone by using tomographic imaging with the data recorded by the pneumatic impact source.
The seismic measurements were conducted using the 3D underground seismic system of the GFZ (Borm et al. 2003) that consists, of a pneumatic impact source, 3-component geophones (28 Hz) for borehole emplacement, and rock anchors. To reduce the effects of the ambient noise, the receivers were placed in the rock anchors at a depth of 2 m (Borm et al. 2003). The receivers were spread over 2 levels in height that were 2 m apart. To further improve the signal-to-noise ratio, the impact source hit 5 times at each source point and the individual hits were stacked. The impact source has a peak force of approx. 4000 N (Borm et al. 2003) and can generate signals with a typical frequency response of 100-2000 Hz in crystalline environments (Richter et al. 2018). For this tomography study, we focused on the partial area with the best ray coverage (

Data selection
The seismic traveltimes were recorded with a 0.25 ms sampling interval and a recording time of 500 ms. For the full data set the normalized mean amplitude spectrum was calculated. We observed a lower frequency range of ∼185-700 Hz within −3 dB (Fig. 2) presumably due to the strong attenuation of clay stone in contrast to crystalline environments (Richter et al. 2018).
With 16 receivers and 19 shot points a total of 304 shotreceiver combinations were acquired. The traveltimes were picked manually for both P-and S-waves. We picked over all three components to ensure to get the right first arrival especially for the S-wave. Fig. 3 shows a receiver gather for receiver 9 in the TT niche. The first arrival of the P-wave (Fig. 3 Fig. 3 Example of a receiver gather at receiver 9 with shot points starting in niche 8 with SP12 to SP31 (except SP23) in gallery 08 ( Fig. 1A along gray dotted line). Red dots indicate traveltime picks for P-waves and blue dots for S-waves. Shown are the three components of the geophone. comp 1: vertical, comp 2: parallel and comp 3: perpendicular to the tunnel wall components, whereas for the S-wave (Fig. 3, blue dot) the quality of the arrival depends on the propagation direction of the wave. For shot point 12-17, located in niche 8, a clear first arrival is seen on component two and for shot point 24-31, located in gallery 08, the first arrival is clearly identified on three.

, red dot) is identical on all
We investigated the data for shear wave splitting by rotation of the components. We hereby focused on source-receiver combinations with ray paths approximately parallel the bedding. A comparison of rotated and unrotated gathers showed that, if quasi-shear waves could be identified, these were the fast phases and consistently polarized parallel to the bedding. We therefore picked the shear-waves on all traces, assuming them to be the fast shear waves when running along the bedding direction.
For the tomography we only use the lower receivers, so that all data refer to a horizontal plane defined by the shot and receiver positions. We also excluded sourcereceiver pairs with small offsets and with ray paths running along on gallery 08 or niche TT/8. For these ray paths, we observed anomalously high apparent velocities, as along these ray paths, waves travel predominantly within the shotcrete at high P-wave velocity of ∼3600 m/s (Guevremont et al. 2000) which is much faster than the adjointing EDZ. Out of 152 possible P-and S-wave picks, 98 (P) and 83 (S) remain for further analyses.

Anisotropic characterization
The anisotropy of Opalinus Clay in Mont Terri is well studied. It is generally described by an oblate ellipsoid with the fast semi-axis parallel and the slow semi-axis perpendicular to the bedding. Our experimental setup is oriented in the horizontal plane of the laboratory and is ∼ 45 • divergent to the symmetry axis of the ellipsoid. To characterize the anisotropy in this section of the ellipsoid, we used the traveltimes of the whole survey and calculated their apparent velocities for straight rays. By fitting an ellipse with a least squares method (Gal 2003 modified) to the azimuthally distributed velocities, anisotropic parameters were estimated, i.e. the fast ( v max ) and slow ( v min ) velocity and the angle of the fast velocity axis. The derived values were then used to calculate the average seismic P-and S-wave velocity anisotropy AvP and AvS: For a better interpretation of the azimuthal distribution of the seismic P-wave and S-wave velocities (Fig. 4), the data were sub-divided depending on their ray paths. The subset 'shaly' included waves that traveled only through the shaly facies (gray dots for source-receiver combinations within the tomography area and black dots for source-receiver combinations from and to the upper part of the pilot survey, Fig. 4) and the subset 'crs' included waves that traveled only through the carbonate-rich sandy facies (gray triangles, Fig. 4). The P-wave velocities of the two subsets 'shaly' (gray/black dots) and 'crs' ( gray triangles) are aligned along an ellipse and the velocities (1) (v max − v min )/v mean . Fig. 4 Scatter plot of the anisotropic characterization in the horizontal plane starting from north ( 0 • ). The apparent P-wave (A) and S-wave (B) velocities were grouped depending on their ray paths: gray triangles-propagation within the carbonate-rich sandy facies and black/lightgray dots-propagation within the shaly facies. The black dots indicate ray paths from and to the upper part of the pilot survey (north of the main fault) through the shaly facies. They are used to increase the azimuthal coverage. Shots parallel to the tunnel wall within the EDZ were excluded. The shown anisotropy ellipses were calculated for the shaly (black) and the carbonate-rich sandy facies (gray) are faster in the NE and SW compared to the NW and SE (Fig. 4A). The resulting characteristics (Table 1) indicate that the different facies are clearly distinguishable.

Tomographic inversion
For the first arrival traveltime tomography we used "simulr16", a 3D simultaneous refraction and reflection seismic tomography by Bleibinhaus (2003). The program is based on simulps12 and simulps13q, which are part of the SIMUL family (Thurber 1983;Thurber and Eberhart-Phillips 1999;Evans et al. 1994).
The inversion program simulr16 is chosen, because of the implementation of anisotropic forward modelling. Thereby, an anisotropic correction along the ray paths is computed in every forward calculation at the beginning of each iteration, resulting in a model with isotropic velocities. But the most important feature of the program is, that the anisotropic parameters can be set independently for each single grid node (Bleibinhaus and Gebrande 2006; Bleibinhaus and Hilberg 2012). In this way, geological layers or features with different anisotropic parameters can be taken into account and thereby a greater degree of heterogeneity can be considered. The forward modeling is carried out with the artificial ray tracer (ART) and pseudo bending (PB) by Um and Thurber (1987), because the anisotropic correction is only available when using ART/PB.
Simulr16 uses a non-linear inversion approach with a damped least square method. The damping should ensure that model changes do not transfer in regions with low ray coverage (Bleibinhaus 2003). The contribution of each grid node to the inversion is decided by the derivative weight sum (DWS), which describes the relative ray density around a grid node (Husen et al. 2000). Thus only grid nodes with a certain threshold will have an impact on the result of the inversion. Because of the non-linear approach ray paths are calculated for every iteration. Together with this calculation the anisotropic correction is applied on the new ray paths. Mind that the inversion does not invert for anisotropic parameters. They are all fix for each inversion node and have to be assigned at the start of the inversion. Therefore three values are needed: azimuth of the symmetry axis, amplitude of anisotropy and dipping angle. To apply the anisotropic adjustment simulr16 iterates over all ray segments and determines the angle between ray segment and the symmetry axis of the anisotropy to identify the value of anisotropy for this ray segment. We modified the code to an approach which leads to the calculation of an isotropic velocity model (mean value of the anisotropic ellipse). The basis for our modification are simplified approximations for weak anisotropy (Thomsen 1986), which are used by Maurer et al. (2006) in a similar approach for moderate anisotropy on wood samples: : angle between symmetry axis and ray segment. v(α) : velocity depending on alpha. v max : velocity of the fast direction. : amplitude of anisotropy.
For numerical reasons we used cosine instead of sine. Furthermore we changed the dependency from v max to the mean velocity v mean of the ellipse, which results in an additional factor of 0.5:

Parameterization
The grid size for the traveltime tomography was defined by a spacing of 2 m in x-, y-direction. For the z-direction a spacing of 4 m is used to include all receiver and shot points in one layer in the vertical direction. In the internal coordinate system the bounding nodes are set to 0 m and 40 m in x-/y-direction and −8 m to 0 m in z-direction. The bounding nodes are linked to adjacent nodes and do not contribute to the inversion. The resulting shape of the grid contains 20 × 20 × 3 nodes. Because of the fixed boundary nodes, only 1 node is calculated in z-direction resulting in a quasi 2D inversion. To decide on a damping value for the inversion, we evaluated multiple single-iteration inversions with different damping values for each iteration. The root mean square (RMS) of the traveltime residuals as an indicator for data variance and the standard derivation of the velocity model as model variance were taken and plotted as trade off curve (Fig. 5). The value with the best improvement of the traveltime residuals and the lowest increase in model variance at the same time is taken for inversion. We decided on six for the P-wave inversion and 35 for the S-wave inversion.

Table 1 Anisotropic characteristics of the shaly and carbonate-rich sandy (crs) facies
Velocities and angle were estimated by fitting an ellipse to the P-and S-wave data (Fig. 4), whereas the seismic velocity anisotropies AvP and AvS were determined according to Eq. (1). The orientation is related to the bedding strike

P-wave tomography
For the parameterization of simulr16, we did not use residual or distance weights, because of the small scale of the experiment. We decided on a value of 10 as a threshold for the derivative weight sum (DWS), which controls which node is inverted. A value of 10 as minimum is needed because of the low ray density and it is used to mask the results of the velocity models. For the significance test we used the 95% f-test. The P-wave tomography was computed with a homogeneous initial velocity model with v = 3.1 km/s which is the weighted mean velocity between shaly and carbonate-rich sandy facies obtained from the anisotropic characterization.
To calculate the isotropic model (mean velocity model) with higher accuracy the anisotropic parameters (AvP in Table 1) were set for every grid node depending on the corresponding geological unit. The velocity of the tomographic result (Fig. 6A) ranges from 2.4 km/s to 3.9 km/s with a mean value of 2.9 km/s for the shaly and 3.5 km/s for the carbonate-rich sandy facies. The DWS (Fig. 6B) shows a higher density of rays near the tunnel walls and at the centre around x = 330 m and y = 580 m, between the diagonal rays in the shaly facies and rays who cross through the carbonate-rich sandy facies, a very low density is observed.

S-wave tomography and v P /v S ratio
The computation of the shear wave tomography resembles that of the P-wave tomography. For the initial model a weighted mean velocity of v = 1.4 km/s was used. The main feature, the transition between the shaly and carbonate-rich sandy facies is recovered as well (Fig. 7A). The mean velocity of the shaly facies is 1.3 km/s, while the carbonate-rich sandy facies has a mean velocity of 1.6 km/s. The ray density is much lower as in the P-wave tomography, which leads to a big hole in the centre of the velocity model. Since both tomographies (P and S) are computed on the same grid, the resulting v P and v S models are used to calculate the v P /v S ratio (Fig. 8) in the tomography area. The v P /v S model is masked with the combined DWS cut off of both inversions. The shaly and carbonate-rich sandy facies can be clearly discriminated by different v P /v S ratios. The v P /v S ratio of the shaly facies is significantly larger (ca. 2.2-2.4) than the v P /v S ratio of the carbonate-rich sandy facies (ca. 1.9-2.0). But higher values contaminate the carbonate-rich sandy facies around x = 325 m and y = 575 m.

Discussion
The anisotropic characterization of the Opalinus Clay shows a good azimuthal coverage for the shaly facies ( > 90 • ) and an underdetermined coverage for the carbonate-rich sandy facies ( ∼ 60 • ) with few data with a strong focus along the fast direction. This limitation affects also the transmission experiment by Schuster et al. (2017). In that experiment, the azimuthal coverage was almost 90 • for the shaly facies, but only ∼ 30 • for the carbonate-rich sandy facies, which only covers for the fast direction. Comparing our results with other studies, which were all carried out within a plane along the anisotropic symmetry axis (Table 2), we can see similar velocities parallel to the bedding. Even perpendicular to the bedding we see a similarity in the velocity although our horizontal plane deviates ∼45 • from the symmetry axis. This is due to the heterogeneities within the shaly facies. Nevertheless our results for v ⟂ P is the highest as expected. Zappone et al. (2021) investigated the shaly facies and main fault right below our own experiment by means of crosshole seismic tomography. Their data lies within the anisotropic symmetry axis and intersects with our data along the 90 • inclination. The v P value of Zappone et al.  (2021) is at the intersection ∼2.5 km/s, what matches our observation (2.551 km/s). The v P /v S ratios after Tatham and Krug (1985) show a clear contrast between sandstone (1.7), carbonates (1.85) and shale (2.2). In this study, the v P /v S ratio of the shaly facies is ∼2.3 and ∼ 2 for the carbonate-rich sandy facies. This seems reasonable, because our geological units are named after the dominant mineral components, but they are still mixed. Therefore, the v P /v S ratio of the carbonate-rich sandy facies is lower than the shaly facies, because its clay-mineral content is low, with a large amount of sand and carbonate.

Conclusions and outlook
The partial data set from the pilot survey, which was acquired for URL scale characterization, was used for tomographic imaging of shaly and carbonate-rich sandy facies of Opalinus Clay. The sparse data acquisition still allowed tomographic inversion and anisotropic characterization. The anisotropic effects on traveltimes were removed within the inversion code simulr16 to a certain degree for the computing of an isotropic velocity model. The high velocity contrast between shaly and carbonate-rich sandy facies can be resolved by tomographic inversion from a homogeneous initial velocity model. Future studies will include a repeat survey with a denser acquisition to improve the azimuthal distribution and the anisotropic characterization. The aspect of monitoring will be analysed in order to investigate potential timelapse effects in the Opalinus Clay caused by operations in the URL (e.g construction of new gallery, injection experiments).