Characterization of heterogeneities in the sandy facies of the Opalinus Clay (Mont Terri underground rock laboratory, Switzerland)

SUMMARY

This contribution is presenting a multidisciplinary investigation of heterogeneities in a clay rock formation, based on seismic tomography, logging and core analysis, as a reconnaissance study for a diffusion experiment. Diffusion experiments in clay rock formations provide crucial experimental data on diffusive transport of radionuclides (RN) in extremely low hydraulic conductivity media. Previous diffusion experiments, conducted, for example, in the Mont Terri underground rock laboratory within the relatively homogeneous shaly facies of Opalinus Clay, and modelling studies of these experiments have demonstrated that the clay rock could sufficiently well be described as a homogeneous anisotropic medium. For other lithofacies, characterized by larger heterogeneity, such simplification may be unsuitable, and the description of heterogeneity over a range of scales will be important. The sandy facies of the Opalinus Clay exhibits a significantly more pronounced heterogeneity compared to the shaly facies, and a combined characterization and RN diffusion study has been initiated to investigate various approaches of heterogeneity characterization and subsequent diffusion in a heterogeneous environment. As an initial step, two inclined exploratory boreholes have been drilled to access the margins of the experiment location. These boreholes have been used to acquire a cross-hole tomographic seismic data set. Optical, natural gamma and backscattering logging were applied and rock cores were analysed. The integrated results of these investigations allowed the identification of an anomalous brighter layer within the investigated area of the sandy facies of approximately 1 m thickness and with its upper bound at roughly 10 m depth within the inclined exploratory wells. Mineralogical analyses revealed only slight variations throughout the rock cores and indicated that the anomalous layer exhibited a slightly higher quartz content, and locally significantly higher calcite contents, accompanied by a lower content of clay minerals. The anomalous layer was characterized by reduced natural gamma emissions, due to the lower clay content, and increased neutron backscattering likely indicating an increased porosity. Seismic P-wave velocities, derived from anisotropic tomography, exhibited a maximal gradient near the top of this layer. The transition from the overlaying darker rock matrix into this layer has been identified as an appropriate location for the setup of a tracer diffusion experiment in a heterogeneous environment.

INTRODUCTION

Numerous countries are considering clay rock formations as potential host rocks for high-level nuclear waste repositories (e.g. Germany, France, Switzerland and Belgium). Recently, the Jurassic Opalinus Clay formation was suggested by the Swiss National Cooperative for the Disposal of Radioactive Waste (NAGRA) as the host rock for the construction of the Swiss deep geological repository (NAGRA 2022). The Mont Terri underground rock laboratory (URL) near St. Ursanne in Switzerland provides the unique opportunity to get direct access to this clay rock formation and to perform experiments and measurements on the rock under in situ conditions (Bossart et al. 2017).

The Opalinus Clay is subdivided into three main lithofacies types, termed shaly, sandy and carbonate-rich sandy facies, differing in the sedimentation environment, accompanying fossils, the mineralogical composition and internal homogeneity (Hostettler et al. 2017). The internal heterogeneities within the main lithofacies types have been the basis for further division into different subtypes, representing, for example, different clay and quartz contents (Lauper et al. 2018). Because hydraulic conductivity is extremely low in Opalinus Clay (10⁻¹²–10⁻¹⁴ m s⁻¹, Houben et al. 2013), molecular diffusion is the main transport process for radionuclides (RN) migrating through the formation. Recently, Lauper et al. (2018, 2021) introduced five subfacies types (SF1–SF5) to allow for a consistent lithological classification of the Opalinus Clay throughout northern Switzerland.

Previous diffusion experiments conducted in Mont Terri (e.g. DI-A, DI-B, DR, DR-B experiments, Van Loon et al. 2004; Wersin et al. 2004; Soler et al. 2008; Gimmi et al. 2014) focused on the comparatively homogeneous and clay-mineral-rich shaly facies. Numerical simulations of tracer diffusion, as observed in these large-scale experiments, demonstrated that for non- and weakly sorbing tracers the treatment of the clay rock as a homogeneous anisotropic porous medium described the diffusion behaviour in the rock sufficiently well (Wersin et al. 2004; Soler et al. 2008, 2019). However, investigations of tracer distribution around the injection intervals (Gimmi et al. 2014) and in cross-hole diffusion studies (Jaquenoud et al. 2021) point towards small but distinct effects of rock heterogeneity on tracer diffusion. This comprises tracer enrichment in certain spots or tracer exclusion from certain layers along the bedding.

Compared to the shaly facies of Opalinus Clay in Mont Terri, clay rock formations considered in other countries as potential candidate sites for nuclear waste repositories may exhibit a more pronounced heterogeneity (e.g. Lower Cretaceous claystones in Northern Germany; Kneuker et al. 2020). A proper conceptualization and numerical treatment of heterogeneity on various scales (nano-/microscopic scale to metre/decametre scale) and upscaling strategies to approximate the effect of small-scale heterogeneities in large-scale models on radionuclide diffusion are therefore an important prerequisite for reliable safety assessment models (Yuan et al. 2022).

Seismic characterization has been part of many experiments in rock laboratories, covering different facies types and a range of scales. On the one hand, the seismic measurements, when integrated in small-scale mechanical investigations of rock samples, enable a more detailed interpretation of spatially sparse information gained from rock cores. Their mechanical properties, as analysed by triaxial deformation tests (e.g. Schuster et al. 2021), can show significant scatter even within single rock formations. On the other hand, large-scale non-destructive geophysical measurements allow the observation of geological processes related to hydrological experiments in heterogeneous rock formations. As an example, seismic tomography imaging is part of the monitoring concept in an experiment, where CO₂ is injected into a fault zone in Opalinus Clay (Zappone et al. 2021). Fault and surrounding rock reactions on leakage are investigated in the context of CO₂ storage. Tomographic inversion under anisotropic conditions is typically applied after performing anisotropic correction of traveltime data, based on the analysis of anisotropy using the azimuthal dependence of average seismic velocities between source and receiver locations (Esefelder et al. 2021; Zappone et al. 2021). The Tournemire Experimental Platform (TEP) in Aveyron, France is in clay rock as well, and is used for experiments related to the assessment of processes in nuclear waste repositories (Boisson et al. 2001). The Toarcian argillite formation surrounding the TEP is characterized by a relatively strong anisotropy of 30 per cent with an almost vertical symmetry axis (vertically transverse isotropy, VTI). Seismic tomography data were acquired between the main tunnel and galleries in a horizontal plane. Here, the inversion of the seismic data could be performed under isotropic conditions (Bretaudeau et al. 2014).

Wawerzinek et al. (2022) analysed the traveltimes of P and S waves from seismic investigations covering the shaly and the sandy facies of Opalinus Clay in the Mont Terri URL. They found pronounced velocity anisotropy (A_v = (v_max—v_min)/v_min) for P- and S-wave velocities (V_P and V_S), particularly for the shaly facies. According to these observations, P-wave anisotropy of 0.261 for the shaly facies contrasts with a significantly lower P-wave anisotropy of 0.094 for the sandy facies. Both values are probably underestimating the ‘real’ anisotropy of the respective facies types as the seismic measurements were performed in a horizontal plane. Due to the inclination of the bedding of the Opalinus Clay at the experiment location, v_min was observed at 45° deviation from the symmetry axis. The results were compared to previous studies and were found consistent considering realistic error margins. While the shaly facies of Opalinus Clay has been investigated in many experiments, less data are available for the sandy facies (Siegesmund et al. 2014; Schuster et al. 2017; Lozovyi & Bauer 2018). Schuster et al. (2017) performed seismic transmission experiments in the sandy, shaly and carbonate-rich sandy facies. A seismic source was activated in Gallery 98 (Ga98, Fig. 1), receivers placed in the Safety Gallery. For the sandy facies, an average P-wave velocity along bedding strike of 3.41 km s⁻¹ was observed. The average P-wave velocities were significantly reduced in the same rock segment for source receiver paths running up to 45° relative to the bedding strike direction. Due to limited azimuthal aperture, anisotropy was not further quantified. High-resolution interval velocity measurements (IVM) along a borehole through the three different facies types of Opalinus Clay revealed very high spatial heterogeneity of apparent V_P in the sandy facies ranging from ≈ 2.5 to ≈ 4.5 km s⁻¹ in the dm to m-scale. The mean values for segments of 1–5 m ranged from 3 to 3.7 km s⁻¹ for V_P perpendicular to the bedding direction (Schuster et al. 2017). Siegesmund et al. (2014) performed seismic velocity measurements at small-scale samples (cubes of 4.3 cm edge length) and estimated average velocity anisotropy (A_vavg = (v_max—v_min)/v_mean) for the shaly facies of 0.25 and for the sandy facies of 0.14–0.16, which translates to slightly higher velocity anisotropy with v_min as a reference. The sandy facies is characterized not only by very high spatial heterogeneity of seismic velocities. The heterogeneity affects also the anisotropy of velocities, which, however, is difficult to quantify in detail due to typically restricted spatial and azimuthal coverage.

The multidisciplinary investigations on the sandy facies of Opalinus Clay in Mont Terri presented here are performed in preparation of a new diffusion experiment in the Mont Terri URL. The so-called DR-D experiment (‘Heterogeneity of sandy facies by geophysical characterization and diffusion studies’) specifically targets a heterogeneous section of the sandy facies to study the impact of rock heterogeneity on RN diffusion. The location of the experiment within the URL is shown in Fig. 1. On the one hand, the results will serve as an initial step for future exploration, and on the other hand, as a proxy for the expectable diffusion behaviour in potentially relevant heterogeneous clay rock formations in countries where clay rock is considered as host rock for nuclear waste disposal. The DR-D experiment will characterize the heterogeneity on various scales. As a first step in this experiment, a rock zone in Gallery 18 in Mont Terri was accessed by two boreholes, BDR-D1 and BDR-D2, and a high-resolution seismic tomography survey was conducted. The seismic investigations represent the largest of the investigated scales (decimetres—decametres). Seismic results were correlated with borehole logging data and information from drill core analyses. The aim of the investigation was to test whether heterogeneities in seismic velocities, detectable by high-resolution seismic tomography, may represent rock sections, which as well exhibit a complex diffusion behaviour. If this approach succeeds, non-destructive seismic tomography may characterize relevant rock formations, and the obtained information may provide initial hints about the expected diffusion behaviour.

SEISMIC TOMOGRAPHY

Acquisition

An overview of the location of the DR-D experiment within the Mont Terri URL is given in Fig. 1. The experiment is situated near the southeastern end of Gallery 18 (Fig. 1a). Two boreholes, BDR-D1 and BDR-D2, were drilled obliquely downward from the floor of Gallery 18. They span an experimental plane which is approximately perpendicular to the bedding of the Opalinus Clay formation (Fig. 1b). The remaining angle between the observation plane and the bedding normal direction is approximately 20°. The RN diffusion setup will be deployed in a third borehole to be drilled between the two exploratory wells. The boreholes were not drilled precisely perpendicular to the bedding plane in order to avoid stability issues during drilling and core retrieval. This will be especially important for the central diffusion borehole. In the following, we will use the term ‘vertical’ indicating the vertical direction within the rotated experimental plane, which is perpendicular to the intersection between experimental plane and gallery surface and almost perpendicular to the bedding planes. The two boreholes are 15.5 m deep and have a diameter of 131 mm. The cross-hole tomographic layout is completed by a receiver profile on the floor of the Gallery 18, between the well heads of both boreholes.

The detailed distribution of seismic source and receiver point locations within the experimental plane is shown in Fig. 2. For the cross-hole part, the impulsive borehole source was operated in BDR-D2 at 30 source point locations. The shots from BDR-D2 were acquired using a 7-level borehole receiver tool in BDR-D1 and 40 surface receivers installed on a profile on the gallery's surface between the two boreholes. During the acquisition, the source was activated on all 30 source point locations. Then, the borehole receiver tool was moved to the next level and the source was activated on all source point locations again. This procedure was repeated as long as needed to acquire seismic data at 15–20 cm receiver intervals along BDR-D1. For the symmetric offset vertical seismic profiling acquisition, the 40 surface receivers acquired all source activations in BDR-D2 during the cross-hole survey. Then, the seismic source was deployed in BDR-D1 and activated at 30 source point locations in that borehole. The resulting data set comprises 599 shot gathers with 47 traces (source in BDR-D2) and 30 shot gathers with 40 traces (source in BDR-D1), respectively. As the fixed surface receivers acquired every source point from the cross-hole survey 20 times (19 times for source point (SP) 30), many traces refer to the same source–receiver combinations.

Data

Fig. 3 shows two exemplary shot gathers of the acquired data after sorting the borehole receiver traces according to the depth of the receivers in the borehole and with one shot recording from the surface profile. The shot gathers show exemplarily the different data quality of the cross-hole data (traces 1–133 for source point 30, and traces 1–140 for source point 46) and the borehole-surface data (traces 134–173 for source point 30, and traces 141–180 for source point 46). The cross-hole data display a lower signal-to-noise ratio, especially at the deepest and shallowest receiver positions, than the borehole-surface data. In addition, the cross-hole data show polarity reversals. These polarity reversals stem from the varying orientation of the seismic borehole source. As an example, traces 50–56 from source gather 46 (Fig. 3) show clearly inverted polarity compared to the adjacent traces. The cross-hole data consist of 20 single gathers of seven traces sorted into one gather, but recorded by 20 different source activations while the source was moving up and down in the borehole. During this mode of operation, the horizontal orientation of the source was not controlled, resulting in these polarity changes along the gathers. The first break arrival times were picked on the data after geometry attribution and sorting into common source point gathers. Thus, the pick data set consists of 60 shot gathers with a total of 2610 traces.

An initial quality control of the traveltime picks is shown in Fig. 4(a). With source–receiver distances larger than 3 m, the average P-wave velocities, when assuming straight rays, range from approximately 2.2 to 3.8 km s⁻¹ with the highest values around 12 m offset. At smaller offsets (below 3 m), lower average velocities were observed, ranging from 1.2 to 3.3 km s⁻¹. The lowest velocities have been observed for source–receiver paths within the excavation damage zone (EDZ) of the gallery. Fig. 4(b) shows the same average P-wave velocities, depending on the azimuth of the straight ray path, with 0° being a path from the source location vertically upward within the observation plane. The slowest average velocities are observed for ray paths propagating vertically within the experimental plane, the fastest average velocities can be attributed to ray paths in horizontal direction parallel to the bedding plane. The azimuthal dependence of average P-wave velocities shows, with some scatter, a clear elliptical form. The Thomsen notation for weak anisotropy (Thomsen 1986) was applied to numerically analyse the azimuthal dependence of average velocities:

V_P(0) is the P-wave velocity perpendicular to the bedding (⁠|$\theta = 0^\circ $|⁠, ‘slow direction’), ε is the P-wave anisotropy and δ is the deviation from an elliptical wave front (Mavko et al. 1998). A minimum root-mean-square (RMS) traveltime fitting for the average values of V_P(0), δ and ε was performed, using a grid-search approach. The resulting best-fitting ellipse is shown in Fig. 4(b) (red line). In order to quantify the uncertainties of the fitting results, the parameter ranges for RMS traveltime errors up to 10 per cent larger than the best fit were computed. Fig. 4(c) shows the best-fitting Thomsen parameters with error bars for a 10 per cent error margin. The best-fitting Thomsen parameters are V_P(0) = 2.63 (± 0.12) km s⁻¹, |$\delta = - 0.01$|(± 0.18) and |$\varepsilon = 0.24$| (± 0.05). This corresponds to a P-wave anisotropy of 24 per cent, translating to approximately 22 per cent average anisotropy, which is slightly larger than previously observed anisotropy on rock samples of the sandy facies (14 per cent–16 per cent according to Siegesmund et al. 2014). The error bars in Fig. 4(c) show that the parameter ε is relatively well constrained, compared to the parameter δ. Comparable observations were made by Le Gonidec et al. (2012) by analysing P-wave traveltime measurements in Opalinus Clay during the excavation of the Gallery 08. The error bar of the velocity V_P(0) is a result of all uncertainties typical for traveltime tomographic investigations (pick uncertainties, heterogeneity of P-wave velocity and anisotropy).

Tomographic inversion of the traveltime data

A Markov chain Monte Carlo (McMC) approach was applied for tomographic inversion of the traveltime data. The approach consists of two main elements: (1) the McMC sampling and (2) the calculation of the (anisotropic) forward problem. A detailed explanation of the McMC scheme is given in Ryberg & Haberland (2018), and the anisotropic extension of the workflow is described in Haberland et al. (2023). Briefly, the model is defined by a set of unstructured points p_i = (x_i, z_i, V_p(0)_i, ε, δ) where x_i and z_i being the spatial coordinates of 2D Voronoi cells, and V_p(0)_i, ε, and δ representing the vertical P-wave velocity, the P-wave anisotropy and the deviation from an elliptical wave front, respectively. The model parametrization is shown exemplarily for two Markov chains in Fig. 5. In order to reduce the number of independent unknown parameters, and considering that the investigation was carried out within one facies type (the upper sandy facies of the Opalinus Clay) characterized by relatively uniform formation history, the parameters ε and δ were assumed to be constant for the entire model. For solving the forward problem, a densely spaced regular 2-D grid was generated and to each gridpoint, the P-wave velocity of the nearest Voronoi cell node was assigned. The same was done for ε and δ, but these consisted of only one Voronoi cell, resulting in homogeneous distributions of ε and δ.

Synthetic traveltimes were computed using a finite-difference (FD) solution of the eikonal equation, extended for anisotropic media by Riedel (2016). The RMS of the difference between the measured and the calculated traveltimes was used as a misfit function in the McMC sampling routine. Thorough and relatively fast sampling of the complete model space was ensured by running 200 Markov chains in parallel. Each chain consisted of 10 000 models, for each Markov chain a random starting model, consisting of a variable number of Voronoi cells with random velocities V_p(0), ε and δ was selected. After the burn-in phase (significant reduction of RMS misfit by convergence to the approximate noise level of the data, 8000 models, see Fig. 5e), the chains were thinned by retaining only every 20th sample of the model space. Finally, 20 000 models, all of them explaining the observed traveltimes equally well, were used for building the final model by locally averaging the velocities from all models, and for performing a statistical analysis to determine the velocity uncertainty/errors (Fig. 6).

Inversion results

The result of the Bayesian traveltime tomography applied to the field data is shown in Fig. 6. The panels in Fig. 6 show V_P(0) (A), the uncertainty of V_P(0) (B), the absolute value of the gradient of V_P(0) (C), and the relative probability of a cell boundary (D). The gradient of V_P(0) and the relative cell boundary probability show almost similar patterns, indicating areas of a high degree of spatially heterogeneous distribution of the P-wave velocity. The inverted velocity model shows the following main features:

1. Lateral and vertical variation of V_P(0), with values ranging between < 2.0 and > 3.0 km s⁻¹, with the highest values in the deep part of the model (below Y = −10 m).

2. Reduced V_P(0) close to the floor of Gallery 18 (within 1 m below the floor), except for a small area in the central area, possibly correlating with the EDZ of the gallery.

3. A moderate high-velocity layer on the left side of the model, at approximately Y = −7 m, intersecting the borehole BDR-D1.

4. Lower V_P(0) close to the boreholes than in the central area of the model, possibly indicating the presence of a damage zone around the boreholes.

The velocity uncertainty (Fig. 6b) is the standard deviation of V_P(0) of all single Voronoi cell models averaged to derive the model in Fig. 6(a). The distribution of the velocity uncertainty shows primarily badly resolved areas of the model. These are (as expected) mainly located at the margins of the model, which are not covered by ray paths between the source and receiver positions in the boreholes and on the gallery floor. Within the well-resolved area of the model, in the central area between the exploration boreholes, higher variance values partly coincide with areas of relatively high-velocity gradient or boundary probabilities. Thus, the uncertainty distribution within the well resolved model can also be used as a detection approach for heterogeneity. The distribution of the gradient and the relative boundary probability (RBP, Figs 6c and d), show maximum values at the same positions, which have been identified in the velocity model (features 1 and 2).

Recovery tests from synthetic data sets

The inversion of synthetic traveltime data, generated for heterogeneous anisotropic velocity models and using the same source and receiver distribution as for the field data was applied in order to identify potential limitations in the interpretation of real data inversions. Synthetic data sets were generated using anisotropic velocity models with varying complexity. The synthetic traveltimes were computed for those source–receiver pairs for which traveltime picks of the field data existed. Fig. 7(a) shows a relatively simple (1-D) model, with a laterally homogeneous P-wave velocity (V_P(0)) distribution, but with a varying vertical gradient, including a velocity inversion (low-velocity zone) between 7 and 9 m depth. The traveltimes were modelled using a FD fast marching method using second-order Godunov schemes (Riedel 2016). Random traveltime jitter (noise), similar in magnitude to the recovered noise level of the real data inversion, was added to the modelled traveltimes (RMS amplitude 0.12 ms). The anisotropy was described by the three Thomsen parameters for weak P-wave anisotropy (V_P(0), ε, and δ), with varying V_P(0) and constant values for ε (⁠|$\varepsilon = 0.25$|⁠) and δ (⁠|$\delta = 0.06$|⁠). It was assumed that under the given conditions in the sandy facies of the Opalinus Clay, the parameters ε and δ, describing the azimuthal dependence of P-wave velocity can be assumed to be constant, reflecting the dominant bedding orientation of the Opalinus Clay, with its rotational symmetry axis almost parallel to the vertical axis of the observational plane. This is the same assumption which is made for the interpretation of the measured field data. Fig. 7(b) shows the result of the inverted velocity distribution (V_P(0)). The parameters ε and δ were inverted for using a single cell model. The inverted values are |$\varepsilon = 0.253$| (± 0.008) and |$\delta = 0.028$| (± 0.023). The first-order ellipticity was almost exactly reproduced by the tomographic inversion, with an error of 1.2 per cent, whereas the deviation from the ellipticity (δ) was not as well reproduced, with an error of 53.3 per cent. These observations are coherent with the error estimations for the range of Thomsen parameters using the straight ray velocity analysis discussed above in the chapter ‘Data’ and shown in Fig. 4.

The difference between synthetic and recovered models of the inverted P-wave velocity distribution is illustrated in Fig. 7(c). The central area of the model between the exploration wells and above 9 m depth was well reproduced with very low-velocity errors of less than 0.1 km s⁻¹. Larger errors are visible at the uppermost part of the area covered by the measurements, right below the receivers on the gallery floor. Here, the simulated EDZ of the gallery was slightly overestimated resulting in too low values of V_P(0), and the low-velocity zone between 7 and 9 m depth was not well resolved, resulting in too high values of V_P(0). The relatively sharp velocity contrast at 9.5 m depth, however, was generally reproduced, with some lateral variability. The parameters variance of the P-wave velocity, velocity gradient and boundary probability, shown in Figs 7(d)–(f), show coherently high values along pronounced heterogeneities within the resolved model area. The variance of P-wave velocity identifies also the areas which are unresolved (or only poorly resolved) at the boundaries of the model and outside the acquisition layout.

Building on the laterally homogeneous 1-D velocity model, a laterally heterogeneous velocity model was generated and the synthetic traveltimes were computed. This model and the inversion results are shown in Fig. 8. The original model was generated using the previously discussed 1-D (laterally homogeneous) velocity model with two types of velocity anomalies added. The first type of anomalies were low-velocity zones around the exploration boreholes, with a radius of 1 m, simulating a hypothetic damage zone of the boreholes. The second type of anomalies was a ‘one-column checkerboard’ pattern in the central area of the model. It consisted of a 2-D sine function describing positive and negative velocity anomalies with a wavelength of 12 m. This results in a high-velocity anomaly around 8 m depth and two low-velocity anomalies around 2 and 12 m depth. The parameters ε and δ were 0.25 and 0.06, respectively, as already discussed in the simpler model above.

The inverted velocity model is shown in Fig. 8(b), and the error of inverted P-wave velocity versus the model is shown in Fig. 8(c). In the velocity model, the upper low-velocity zone and the high-velocity anomaly in the central part of the model could be identified, but with slightly different shapes. The sharp upper boundary of the low-velocity anomaly in the synthetic model was not reproduced by the tomography, and also the layer boundary right below the high-velocity anomaly (at 9.5 m depth) was not well detected by the tomography. The low-velocity anomalies around the boreholes have been detected in the upper part of the model (down to about 6 m depth). Below 6 m depth, only some faint and blurred velocity reduction around the boreholes may be interpreted as a result of the low-velocity zones around the boreholes in the forward model. Especially the lower part of the model (below Y = −9 m) close to the left source and receiver borehole was not well resolved. The large velocity gradient at around Y = −9.5 m in the model was not reproduced by the tomographic inversion close to the left borehole. Close to the right borehole, the same feature is visible, however with some dislocation of 0.5–1 m. The distribution of velocity errors (Fig. 8c) shows that in the central area of the model, the inversion result has mostly underestimated the ‘real’ P-wave velocities, whereas along the boreholes, velocities were overestimated. The parameters ε and δ, again, have been inverted using one single cell for the whole model as for the previously discussed model. The inverted values are |$\varepsilon = 0.271$| (± 0.011) and |$\delta = 0.068$| (± 0.030).

As for the first model, the standard deviation of V_P(0), the gradient of V_P(0), and the RBPs show coherently areas of enhanced heterogeneity within the well-resolved part of the model. The boundaries of the resolved high- and low-velocity anomalies were detected by high gradient and boundary probability values. A high V_P(0) standard deviation was also imaged along the boreholes, indicating the high uncertainty of inverted velocities here due to insufficient resolution of the damage zones around the boreholes.

BOREHOLE LOGGING AND CORE CHARACTERIZATION

Borehole logging

Subsequently after drilling the boreholes BDR-D1 and BDR-D2, several borehole loggings were acquired. An optical logging was used for a structural analysis of the rock. The Schmidt plot confirmed the inclination of the clay rock bedding of about 55°, dipping towards 145° (SE-SSE). Otherwise no significant structural features were detected.

A spectral gamma log revealed no structural correlations between the two boreholes or other obvious rock features. The overall activity was not sufficient for isotope specific signals beyond noise level. However, the overall gamma counts indicate some systematic variation (Fig. 9, left-hand graph, API).

Radioactive loggings of the two boreholes with a ¹³⁷Cs source for backscattered gamma detection and a ²⁵²Cf source for backscattered neutron–neutron detection were performed in triplicate. The backscattered gamma counts were recalculated to a rock density (g cm⁻³, GD, Fig. 9). The density values showed a slightly increased scattering in the EDZ of Gallery 18. Below that, they did not indicate any significant or systematic density variations. Neutron–neutron backscattering counts are depicted in Fig. 9 as well (NN). Interestingly, the backscattered neutron counts are perfectly anticorrelated with the natural gamma counts, although the signal is noisy. As neutron backscattering is correlated with the abundance of light elements (especially H), the neutron backscattering signal may in some cases directly be correlated with pore water content. However, in clay mineral-rich rocks like the Opalinus Clay, the OH-groups of the clay minerals contribute to neutron–neutron backscattering and may even outcompete the pore water contribution.

Drill-core analyses

The extracted drill-cores were polished, photographed and visually classified into the following categories: dark-grey clay matrix, light-grey layers and concretions, turbated light-grey matrix and brownish concretions. This was done by assessing the dominant phases in each square of a rectangular grid with net width of 2·2 cm². Multiple repetitions of this subjective interpretation revealed an uncertainty of ± 10 per cent for this assessment. Results of this simple categorization are shown in Fig. 9 and used for a comparison between the two drill cores in Fig. 10. On the one hand side, they demonstrate the variability of the rock perpendicular to the bedding. On the other hand side, they highlight the homogeneity of the rock along the bedding planes on a 10–20 m scale, as demonstrated with the strong correlation between the rock categories identified in the two drill cores from 5.0 to 15.5 m. Even in the small-scale image of the two drill cores depicted in Fig. 10, specific brighter layers can be seen to occur in both cores at equal depth. An especially prominent feature, a fossil-rich layer at 7.25 m is marked in red in Fig. 10. A series of drill-core images of BDR-D2 from 8.0 to 12.5 m is shown in Appendix  A.

In the direction perpendicular to the bedding, the rock shows a pronounced heterogeneity. The majority of the drill-cores appears as ‘homogeneously heterogeneous’, with a dominating dark-grey clay matrix, which is intercalated by 10–50 per cent light-grey or whitish layers and lenticular structures, varying in thickness, size and shape, from 1 mm thin layers to rounded structures with several cm in diameter (cf. images of core BDR-D2 in Appendix  A). Some of the features found in the BDR-D drill-cores, especially the ‘brownish concretions’, can even be traced in niches 3 and 4 of Gallery 18, approximately 40 m to the southwest of the DR-D experiment (see Fig. 1b for location). A similar fossil-rich layer as was detected at 7.25 m, is indicated as coquinas, and identified as pelecypods, by Hostettler et al. (2017) from the BDB-1 core, approximately 8 m below the transition from the Passwang formation to the Opalinus clay. This is suggesting that this layer may even have a lateral extent of > 150 m.

A layer that is clearly distinct from the general texture was observed between 10.1 and 11.1 m well depth (WD). Here, the rounded light-grey or whitish textures are dominating, and the usual bedding direction is hard to trace. This disordered light-grey or whitish layer correlates perfectly with the prominent signatures observed in the natural gamma and the neutron backscattering logs at the same depth (cf. API and NN signals in Fig. 9). Relating this observation to the description of the sandy facies by Hostettler et al. (2017), the light-grey layer may likely be associated with the concretionary limestones described there, which are characterized as sandy limestones in a matrix of sandy or silty clay stone. With a purely visual assessment, the bright region from 10.1 to 11.1 m can be classified to the subfacies SF3, according to the nomenclature scheme of Lauper et al. (2018, 2021). The rest of the bore cores optically shows great concordance with the SF2.

To get a more profound insight into the nature of this layer compared to the general rock matrix, a mineralogical characterization of the BDR-D1 and BDR-D2 cores was performed. To this end, approximately every 25 cm a 2 cm thick slice of the cores was cut out, and divided in two halves. One of these halves was crushed and milled to powder in a ball mill (Fritsch Pulverisette) and investigated by X-ray powder diffraction (XRD) employing a Bruker D8 Advance diffractometer without further pre-treatment. Diffractograms were quantitatively analysed via Rietveld refinement using the Bruker Topas 7.0 software. An exemplary refinement and a description of the procedure are shown in Appendix  B.

XRD revealed that the clay rock consists of quartz, the carbonates calcite, dolomite and siderite, the clay minerals kaolinite, illite and chlorite (trace amount of illite-smectite mixed layers), k-feldspar/orthoclase and albite, with occasional minor amounts of pyrite and apatite, which is in agreement with previous investigations (Bossart et al. 2017). For a more robust classification, categories of minerals were defined and several mineral fractions combined. A ternary diagram shows the mineral distribution of the investigated samples (Fig. 11, left) categorized as quartz, carbonates and clay minerals, normalized to 100 per cent total. The plots do not include the fossil layer, which stands out with an extremely high carbonate content. The approximate mineral compositions assigned to the SF2 and SF3 subfacies, Lauper et al. (2018, 2021) are indicated in green and red, respectively. The visual classification holds true, and the identified mineral compositions coincide with the reported composition ranges. The brighter samples from the range of 10.0–11.0 m WD show a slightly higher quartz- and carbonates content, and are correspondingly depleted in clay minerals. The samples from the fossil layer and the layer adjacent to it show an above average carbonates content, ranging up to more than 65 per cent in both cores. A more precise mineralogical composition is displayed in Fig. 11 on the right with a side-by-side comparison of the two cores at similar depths. Selected diffractograms are depicted in Appendix  B. In general, the most abundant mineral is quartz (30–43 per cent), followed by kaolinite (13–21 per cent), illite (10–18 per cent) and chlorite (5–11 per cent). Of the carbonates, calcite is most common (4–13 per cent) followed by siderite (3–6 per cent) and dolomite (1–4 per cent). The remaining minerals present are feldspars (8–13 per cent) and trace amounts of pyrite and apatite. An important outcome of the mineralogical characterization is that brighter and darker layers in the clay rock cannot be generally attributed to a certain mineral composition. However, locally the bright layer contains as expected areas which are strongly enriched in carbonate minerals (e.g. samples at 10.25 m, especially D2), and correspondingly depleted in clays and feldspars. This effect is especially pronounced in the fossil layer, where calcite amounts to 65 per cent, while illite contents are as low as 5 per cent. Only in this coquina layer, a relevant pyrite content has been observed with 2 per cent. However, this effect is rather constricted to single lenses or concretions and does not characterisz the whole layer.

The lower clay and feldspar content likely explain the negative peak in the natural gamma emission log (Fig. 9, API). The increased neutron–neutron backscattering (NN) in the bright layer cannot be explained by the mineralogical composition, it is likely due to the presence of pore water, and may thus indicate an increased porosity within this layer. Due to the drilling, transport and storage conditions, drying of the samples could not be prevented prior to measurements of the water content via thermogravimetric analysis. With this analysis, a homogeneous pore water content of 1.0–1.5 per cent has been observed across both cores with no significant deviation in the light-grey area of interest.

INTEGRATION OF RESULTS AND DISCUSSION

Fig. 12 shows the results of seismic tomography between the boreholes and the floor of Gallery 18 together with logging data from the boreholes. The logging data comprise the natural gamma counts, gamma density and neutron–neutron backscattering. As described above, the natural gamma and the neutron–neutron backscattering logs show only little systematic variation along the borehole trajectories with an exception between 9.5 and 11.0 m true vertical depth within the experimental plane (TVD). Lower values of the natural gamma counts correlate with higher values of neutron backscattering. It would be desirable to compare these punctual and high-resolution logging observations with the results of the 2-D seismic tomography, expanding the observations into the space between these wells. However, an analysis of the error of the tomography result of V_P(0) (see Fig. 6b) shows that it is significantly larger very close to the borehole trajectories (source and receiver profiles) than it is in the centre of the investigated area. Therefore, the seismic velocities along the borehole trajectories should be interpreted with care. On the other hand, the natural gamma and neutron gamma logs along both boreholes show similar main features at the same depth levels, suggesting that the relevant properties, causing the observed signatures, are aligned along layers parallel to the bedding strike, which is horizontal within the observation plane.

The P-wave velocities were extracted from tomographic inversion along trajectories, which were shifted parallel from the observation boreholes by 5 m towards the centre of the investigated area for further comparison between seismic velocities and logging results. The shifted trajectories are indicated by thin lines in Fig. 12. In Fig. 12, the depth profiles of P-wave velocity V_P(0) and the absolute values of the velocity gradient along trajectories parallel to BDR-D1 and BDR-D2, respectively, are compared on the left- and right-hand sides of the tomographic image. The focus is on the borehole logging and tomography results below 3 m depth (the approximate penetration depth of the EDZ of Gallery 18 and related scattering of results in borehole logging). From the natural gamma and the neutron backscattering logs, two distinct zones can be identified: zone 1 has high average natural gamma counts (API about 100 counts) and low neutron–neutron backscattering counts (NN about 5 counts) along the larger part of the boreholes. Zone 2 has lower average natural gamma counts (API about 80–90 counts) with higher neutron backscattering counts (NN up to about 6 counts) restricted to a smaller part of the boreholes between 9.5 and 11.0 m depth.

Along a trajectory running parallel to borehole BDR-D1, P-wave velocity is increasing from roughly 2.7 km s⁻¹ at 3 m WD to highest values of V_P(0) (V_P(0) about 2.9 km s⁻¹) within the well-resolved part of the model and the velocity gradient is at its maximum near the top of zone 2. This may be an indication that the different signatures in the borehole logging (and visual core analysis) and the related small variations in mineralogical composition affect the elastic properties of the formation here and thus result in a high-velocity anomaly between 10.5 and 12 m depth. In the tomographic image, below 12 m depth, a weak indication of reduced seismic velocity, compared to the zone above 12 m depth, may be identified. In traveltime-based tomographic inversion, such ‘velocity inversions’ (low-velocity zones between layers of higher seismic velocity) are difficult to detect, which is also demonstrated by the recovery test as shown in Fig. 8. In that test, a velocity inversion at 9 m depth could not be reproduced by tomographic inversion using synthetically simulated traveltime data. It is therefore reasonable to assume that seismic velocities are actually lower at 12 m depth and below than shown in the tomographic result, which would then correlate well with the signatures of the natural gamma and neutron backscattering logs. It is also interesting to note that within around 4 m lateral distance to the borehole BDR-D1, the tomographic result does not show any indication of a transition to high P-wave velocities at the top of zone 2, which is the case for the remaining part of the model. This feature is similar to the unresolved reconstruction of a high-velocity gradient close to the left borehole in the 2-D heterogeneous recovery test shown in Fig. 8(b). In conclusion, it is indeed possible that the tomographic inversion has not been able to detect high seismic velocities within zone 2 close to borehole BDR-D1 due to limited spatial coverage and that the ‘high-velocity layer’ within zone 2 probably extends further towards borehole BDR-D1 than shown in the tomographic result.

In a detailed analysis of a 27.6 m long drill-core covering all three main facies types of Opalinus Clay, including the lower sandy facies (Lauper et al. 2018), high intra-facies heterogeneity was observed, and the lithofacies classification was refined. The lower sandy facies was subdivided into two subfacies types, ‘lenticular’ and ‘sandy’ subfacies, with relatively low (20–40 per cent) and rather high (40–50 per cent) quartz contents, respectively. The sandy subfacies was characterized by higher P-wave velocities and gamma density than the lenticular subfacies. Such high variations in the quartz contents could not be confirmed by our mineralogical analyses within the upper sandy facies, where rather low contrasts in quartz contents but larger variations in calcite contents were found. Taking into account these observations, we attributed the logging and seismic velocity signature within zone 1 to sandy–silty clay-stone (slightly reduced quartz contents), and the signatures of zone 2 were attributed to the concretionary sandy limestone.

CONCLUSIONS AND OUTLOOK

In this study, seismic cross-hole tomography, borehole logging and core analysis results were combined to characterize clay rocks within the upper sandy facies of Opalinus Clay in the Mont Terri URL. The main task of this study, integrating a range of scales and geophysical and structural investigations, was to identify an area of pronounced heterogeneity as a most suitable location of the core activity in the planned DR-D diffusion experiment. The diffusion experiment will investigate the impact of lithofacies variability on the migration of RN in clay rock by injecting a tracer cocktail from a borehole to be drilled into the investigated experimental plane between the two existing exploratory boreholes. The highest degree of variability in lithofacies will, according to the characterization results, be found within the experimental plane at approximately 9.5 m distance from the floor of Gallery 18, between the well heads of boreholes BDR-D1 and BDR-D2. This position, which corresponds to approximately 10.1 m WD (well depth) along the trajectories of the exploratory boreholes, has been identified as the top of an anomalous, approximately 1–1.5 m thick zone with lower natural gamma activity, higher neutron backscattering, higher seismic P-wave velocity, a slightly higher quartz and locally higher carbonate mineral content, and correspondingly lower clay and feldspar contents, compared to the rest of the investigated rock mass. While the logging and structural investigations along the boreholes have provided precise and spatially highly resolved information along the borehole trajectories, the seismic cross-hole tomography has provided a 2-D image of heterogeneous structure between the exploratory wells. Within the range of resolution of the seismic tomography, it confirmed the lateral extension of the anomalies identified by the point measurements within the boreholes. Therefore, it can be concluded, that the deployment of a 1 m injection interval around the upper bound of the anomalous layer will provide ideal experimental conditions for investigating the impact of the presence of various lithofacies subtypes within the Opalinus Clay on radionuclide diffusion.

ACKNOWLEDGEMENTS

This work was funded by the Federal Ministry of Education and Research (grant 02 NUK 053 C) and the Helmholtz Association (grant SO-093). The authors thank the team of swisstopo, in particular David Jaeggi, for their support in all field experimental requests in Mont Terri, Schützeichel GmbH & Co. KG for leading a successful drilling campaign, terratec Geophysical Services GmbH & Co. KG and BBi Brunnen- und Bohrlochinspektion GmbH for the logging of the boreholes as well as Solexperts AG for casing the two boreholes for seismic measurements. We are sincerely grateful for the seismic data acquisition performed by Hendrik Albers, Ulf Nowak, Friedhelm Schulte and Torsten Tietz (BGR). Seismic data visualization and first break picking was done using the Claritas seismic software which is made available to the GFZ via an academic license by Petrosys. The Generic Mapping Tools (GMT) package was used for data and model visualization (Wessel & Smith 1998). The constructive remarks of an anonymous reviewer have helped improve part of the manuscript which we gratefully acknowledge.

DATA AVAILABILITY

The data acquired within this study are available to members of the Mont Terri Consortium.

References

APPENDIX A. IMAGES OF THE POLISHED CORE BDR-D2 FROM 8.0 M DOWN TO 12.6 M (WELL DEPTH)

APPENDIX B. X-RAY DIFFRACTOGRAMS OF DRILL CORE SECTIONS AT VARIOUS DEPTHS AND AN EXEMPLARY RIETVELD REFINEMENT

BDR-D1:

BDR-D2:

Rietveld refinement of core BDR-D2 at 10.0 m WD:

For quantitative analyses of the recorded diffractograms, Rietveld refinements were performed. First, all major components were identified using the program code Diffrac.EVA V5.0 and the crystallographic open database (COD). Next, the CIF files of the identified mineral phases were loaded into the program code TOPAS-7 for the Rietveld refinement. To ensure comparability, a single procedure has been established during preliminary analyses, to fit all samples with a consistent set of adjustable parameters. In this evaluation, the lattice parameter (lengths and angles) of all crystals were allowed to change by 3 per cent and a variable underground is implemented as a background correction. The minimum possible particle size was identified to be the most critical parameter inducing the biggest changes in the quantification. Therefore, each sample was fit in triplicate with three different lower limits for the adjustable particle size of 30, 100 and 200 nm. For the final composition, the average and the standard deviation of these three fits was taken. During this whole procedure, the goodness of the fit, a parameter determining the quality of the fit provided by the TOPAS software, was monitored to maintain adequate quality (<4). Additionally, synthetic mixtures of selected minerals in different compositions were prepared, measured and fit, to confirm the quality of the fitting procedure.

Author notes