Geophysical Journal International

SUMMARYWe present a time-domain distributional finite-difference scheme based on the Lebedev staggered grid for the numerical simulation of wave propagation in acoustic and elastic media. The central aspect of the proposed method is the representation of the stresses and displacements with different sets of B-splines functions organized according to the staggered grid. The distributional finite-difference approach allows domain-decomposition, heterogeneity of the medium, curvilinear mesh, anisotropy, non-conformal interfaces, discontinuous grid, and fluid-solid interfaces. Numerical examples show that the proposed scheme is suitable to model wave propagation through the Earth, where sharp interfaces separate large, relatively homogeneous layers. A few domains or elements are sufficient to represent the Earth’s internal structure without relying on advanced meshing techniques. We compare seismograms obtained with the proposed scheme and the spectral element method, and we show that our approach offers superior accuracy, reduced memory usage, and comparable efficiency.

SUMMARYCorrecting atmospheric effects in measurements of vertical acceleration is challenging, particularly at frequencies above 0.3 mHz. Corrections based on four-dimensional atmosphere models should be the most effective, but are limited in temporal and spatial resolution. So these models are commonly combined with a constant admittance between air pressure and gravity, though the correction can be further improved by allowing a frequency-dependent admittance. This paper studies the effectiveness of such procedures given that air pressure admittance varies with time because of variability in the underlying local atmospheric process, and finds that admittances estimated by cross-spectral analysis gives better results than a constant admittance does. The best results are achieved by combining a 4D model with admittances estimated from time series that cover the event of interest.

SUMMARYDike propagation is a mechanism for more rapid ascent of felsic magmas through the crust than is possible via diapirs or percolative flow. As it ascends, the magma undergoes complex physical and chemical transformations induced by decompression and cooling. These processes dramatically change the magma density and viscosity, which in turn affect magma ascent rate and the depth at which the dike arrests. We present a mathematical model of dike propagation for silicic magmas taking into account the presence of multiple volatile species (H2O and CO2), bubble growth, heat advection and loss, crystallization and latent heat release. We consider conditions for dikes associated with porphyry ore deposits, which may represent an endmember in rapid ascent of felsic magmas from depth. In particular, we simulate the propagation of dikes launched from a deep (900 MPa), volatile-saturated magma source, testing the effects of the magma H2O/CO2 content, temperature, and mass on its ascent rate and final emplacement depth. The model predicts short ascent times (hours to days), with a large increase in viscosity at shallow depth, leading to stagnation and solidification of the dike. Higher initial water content, higher temperature and larger mass of the magma in the dike promote faster propagation and shallower arrest. Volatile loss from ascending magma remains limited until the stagnation depth, providing a potential mechanism for transfer of deep volatiles to hypabyssal blind intrusions associated with porphyry ore deposits. Our findings are applicable to the problem of silicic magma ascent through the crust more generally.

SUMMARYImplementations of Markov chain Monte Carlo (MCMC) methods need to confront two fundamental challenges: accurate representation of prior information and efficient evaluation of likelihood functions. The definition and sampling of the prior distribution can often be facilitated by standard dimensionality-reduction techniques such as Principal Component Analysis (PCA). Additionally, PCA-based decompositions can enable the implementation of accurate surrogate models, for instance, based on polynomial chaos expansion (PCE). However, intricate geological priors with sharp contrasts may demand advanced dimensionality-reduction techniques, such as deep generative models (DGMs). While suitable for prior sampling, these DGMs pose challenges for surrogate modeling. In this contribution, we present a MCMC strategy that combines the high reconstruction performance of a DGM in the form of a variational autoencoder (VAE) with the accuracy of PCA-PCE surrogate modeling. Additionally, we introduce a physics-informed PCA decomposition to improve accuracy and reduce the computational burden associated with surrogate modeling. Our methodology is exemplified in the context of Bayesian ground-penetrating radar (GPR) travel-time tomography using channelized subsurface structures, providing accurate reconstructions and significant speed-ups, particularly when the computation of the full-physics forward model is costly.

SUMMARYThis study focuses on the 3-D velocity structure and thickness of ∼7 Myr-old oceanic crust surrounding borehole 504B, located ∼235 km from the intermediate-spreading Costa Rica Rift (Panama Basin). It investigates how well seismic structure determined by 3-D tomography compares with actual lithology and, consequently, what the origin and cause might be of an amplitude anomaly, the 2A Event, that is observed in multichannel seismic data. Our P-wave model shows an ∼0.3 km-thick sediment layer of velocity between ∼1.6-1.9 km s−1 (gradient 1.0 s−1), bound at its base by a velocity step to 4.8 km s−1 at the top of oceanic crustal Layer 2. Layer 2 itself is subdivided into two main units (2A and 2B) by a vertical velocity gradient change at 4.5 km depth, with a gradient of 1.7 s−1 above (4.8-5.8 km s−1) and 0.7 s−1 below (5.8-6.5 km s−1). The base of Layer 2, in turn, is defined by a change in gradient at 5.6 km depth. Below this, Layer 3 has a velocity range of 6.5-7.5 km s−1 and a gradient of ∼0.3 s−1. Corresponding S-wave igneous layer velocities and gradients are: Layer 2A, 2.4-3.1 km s−1 and 1.0 s−1; Layer 2B, 3.1-3.7 km s−1 and 0.5 s−1; Layer 3, 3.7-4.0 km s−1 and 0.1 s−1. The 3-D tomographic models, coupled with gravity modelling, indicate that the crust is ∼6 km thick throughout the region, with a generally flat-lying Moho. Although the P- and S-wave models are smooth, their velocities and gradients are remarkably consistent with the main lithological layering subdivisions logged within 504B. Thus, using the change in velocity gradient as a proxy, Layer 2 is interpreted as ∼1.8 km thick and Layer 3 as ∼3.8 km thick, with little vertical variation throughout the 3-D volume. However, the strike of lateral gradient variation is not Costa Rica Rift-parallel, but instead follows the orientation of the present-day adjacent Ecuador Rift, suggesting a reorientation of the Costa Rica Rift spreading ridge axis. Having determined its consistency with lithological ground-truth, the resulting P-wave model is used as the basis of finite difference calculation of wave propagation to find the origin of the 2A Event. Our modelling shows that no distinct interface, or transition, is required to generate this event. Instead, it is caused by averaging of heterogeneous physical properties by the seismic wave as it propagates through Layer 2 and is scattered. Thus, we conclude that the 2A Event originates and propagates exclusively in the lower part of Layer 2A, above the mean depth to the top of the dykes of Layer 2B. From our synthetic data we conclude that using the 2A Event on seismic reflection profiles as a proxy to determine the Layer 2A/2B boundary's depth will result in an overestimate of up to several hundred metres, the degree of which being dependent on the specific velocity chosen for normal moveout correction prior to stacking.

SummaryA large non-double-couple component of a tectonic earthquake indicates that its rupture likely was complex and likely involved multiple faults. Detailed source models of such earthquakes can add to our understanding of earthquake source complexity. The 2007 Martinique earthquake in the Caribbean Sea is one of the largest recent earthquakes with a known large non-double-couple component. It was an intermediate depth intraslab earthquake within the South American plate where it is subducting beneath the Caribbean plate. We applied potency density tensor inversion (PDTI) to teleseismic P waves generated by the 2007 Martinique earthquake to model its source processes and focal mechanism distribution. We identified two focal mechanisms: a strike-slip mechanism with a north–south tension axis (T axis), and a down-dip extension (DDE) mechanism with an east–west T axis. Rupture by the DDE mechanism was predominant in the northern part of the source region and strike-slip rupture in the southern part. These two focal mechanisms had approximately parallel pressure axes (P axes) and approximately orthogonal T axes. The seismic moments released by both types of rupture were almost equal. These results indicate that the 2007 Martinique earthquake had a large non-double-couple component. We identified five sub-events with two predominant directions of rupture propagation: two strike-slip sub-events propagated to the southeast and three DDE sub-events propagated to the east. Although the directions of propagation were consistent for each focal mechanism, each sub-event appears to have occurred in isolation. For example, the rupture of one DDE sub-event propagated from the edge of the source region back towards the hypocentre. Complex ruptures that include multiple sub-events may be influenced by high pore fluid pressure associated with slab dehydration. Our results show that PDTI can produce stable estimates of complex seismic source processes and provide useful information about the sources of complex intermediate depth intraslab earthquakes for which fault geometry assumptions are difficult.

SummaryContinuous Plinian eruptions often excite atmospheric modes of ∼3.7 and ∼4.4 mHz, which are observed as harmonic oscillations of ionospheric total electron content (TEC) by global navigation satellite system (GNSS) receivers. Such TEC oscillations started shortly after the great eruption of the Hunga-Tonga Hunga-Ha'apai (HTHH) submarine volcano at ∼4:14 UT, on 2022 Jan. 15. Here I analyze GNSS data at stations within ∼4000 km from the volcano to study temporal and spatial distribution of such atmospheric modes. Strong ∼3.7 mHz TEC oscillations in near fields started shortly after the eruption onset and propagated outward with the sound speed from HTHH. Later such TEC oscillations became strong again with the amplitude peak at the distance ∼1400 km from HTHH. Such far field oscillations occurred also above New Zealand and the Solomon Islands, ∼3000 km from HTHH. Their amplitudes seem correlated with those of the 0S29 solid earth mode, suggesting that vertical surface vibrations underneath may play a role in maintaining the atmospheric mode. Onset of the far field TEC oscillations are synchronized with the local sunrise, possibly controlled by diurnal changes in the ionospheric electron density.

SummaryEarthquake source parameters can be estimated using seismological observations, but the identification of the fault responsible is often complicated by location uncertainties and the inherent ambiguity between nodal planes. Satellite InSAR can be used to observe ground deformation and model fault geometry but is limited by climate conditions (water vapour) and ground coverage (dense vegetation). In the tropics, the atmosphere is dynamic and most regions are densely vegetated, making detecting co-seismic deformation challenging. Here we perform a systematic inspection of co-seismic interferograms from Sentinel-1 SAR images, to assess their suitability for detecting co-seismic deformation in Costa Rica. Using data from the seismological network, we target seven earthquakes between 2016 and 2020 with depths ≤ 20 km and magnitudes Mw5.3 to Mw6.2 . For each event, we use the seismic parameters to compute line-of-sight displacements for ascending and descending geometries and for both nodal planes and generate 12 and 24-day co-seismic interferograms where available. We obtain interferograms with co-seismic displacement signals for three of the seven earthquakes. We invert the geodetic data to retrieve the earthquake source parameters but the lack of offshore geodetic coverage causes tradeoffs between parameters and large uncertainties. The Jacó and Golfito earthquakes likely occurred on the subduction interface and the geodetic locations were 6 − 9 km closer to the coast than previous seismic estimates. The Burica earthquake occurred on a shallow steeply dipping thrust fault in the outer forearc. For the other earthquakes, no co-seismic deformation was detected due to atmospheric noise or poor coherence. These results demonstrate the suitability of 12-day Sentinel-1 interferograms for monitoring shallow earthquakes of magnitude >Mw5.7 in Central America. This approach can be used to begin a surface deformation catalogue for the region, which will ultimately help improve the understanding of active deformation processes and improve hazard maps.

SummaryWe perform seismic and ultrasonic measurements in carbonate and shaley sandstone samples as a function of differential pressure. The velocities show a strong frequency and pressure dependence. The dispersion disappears with increasing pressure and the squirt flow in turn inhibits the pressure dependence. To model these effects, we combine the Gurevich's squirt-flow model with the Mori-Tanaka scheme and the David Zimmerman model, and extend it with third-order elastic constants, to obtain a frequency-dependent acoustoelasticity model. Comparisons between measurements from this study and literature and modeling results show that the P-wave velocity increases nonlinearly first and then nearly linearly, dominated by crack closure and acoustoelasticity, respectively. The pressure dependence of wave velocities is reduced by liquid substitution and further by the squirt-flow mechanism. The effects of fluid properties and crack closure on P-wave velocity decrease with differential pressure. The results will feed a new model and help better understanding the wave propagation in prestressed rocks at different scales.

SummaryDistributed Acoustic Sensing (DAS) is becoming increasingly popular in microseismic monitoring operations. This data acquisition technology converts fiber-optic cables into dense arrays of seismic sensors that can sample the seismic wavefield produced by active or passive sources with a high spatial density, over distances ranging from a few hundred meters to tens of kilometers. However, standard microseismic data analysis procedures have several limitations when dealing with the high spatial (inter-sensor spacing up to sub-meter scale) sampling rates of DAS systems. Here we propose a semblance-based seismic event detection method that fully exploits the high spatial sampling of the DAS data. The detector identifies seismic events by computing waveform coherence of the seismic wavefield along geometrical hyperbolic trajectories for different curvatures and positions of the vertex, which are completely independent from external information (i.e. velocity models). The method detects a seismic event when the coherence values overcome a given threshold and satisfies our clustering criteria. We first validate our method on synthetic data and then apply it to real data from the FORGE geothermal experiment in Utah, USA. Our method detects about two times the number of events obtained with a standard method when applied to 24h of data.

SummaryClairaut's theory that relates the Earth's oblate figure and internal ellipticity to its gravity under rotational-hydrostatic equilibrium has reigned classical geodesy over the centuries. In this paper, we (i) derive from first principles the classical Clairaut's theory for the polar oblateness of a rotating planet under axi-symmetric rotational-hydrostatic equilibrium, and (ii) extend the development to the triaxial case for the equatorial ellipticity of a tidally locked synchronous-rotating moon under rotational-tidal-hydrostatic equilibrium. Typical derivations of the classical Clairaut's theory presented in the literature being rather laborious even to first order, we instead exploit two concise forms of methodology: the gravitational multipole formalism on the physics side, and the Jacobian determinant for the Clairaut coordinate transformation on the mathematics side. The outcome is a logical and straightforward derivation of Clairaut's theory to first order in its entirety, encompassing all the equations and related formulas in geodesy bearing Clairaut's name. That further allows a natural extension to a tidally locked moon. In particular it is demonstrated that the same Clairaut's differential equation applies to both cases governing both the polar oblateness and the equatorial ellipticity.

SummaryFor site characterization, the elliptic particle motion of Rayleigh waves and its frequency dependence is a well-known property that aroused less interest than the frequency dependence of the phase velocity. More than fifty years ago, ellipticity was already recognised as providing information independent from phase velocity, despite the difficulties inherent to its accurate and precise measurement. Several techniques were developed during the last two decades to extract the ellipticity curve from ambient vibration recordings, with a single three-component (3C) station, with pairs of 3C stations and more recently with 3C arrays. The latter has the advantage over the other approaches that the sign of the ellipticity can be retrieved. Moreover, higher order mode separation is possible under certain conditions. Nevertheless, Rayleigh Three-component BeamForming (RTBF) proposed by Wathelet et al. (2018) encounters difficulties in the presence of significant levels of incoherent noise when the true ellipticity is vanishing or when it has a high absolute value. In this work, the analytical expressions of the beam power for a single source wavefield are revised under more realistic assumptions for the incoherent noise azimuthal distribution. The proposed model also includes an asymmetric distribution of the incoherent noise between vertical and horizontal components, which was not the case in the original publication. Switching from ellipticity to angular ellipticity drastically simplifies the formalism. Moreover, it naturally leads to a new steering matrix (All-component ellipticity steering) which solves the limitation around zero and infinity observed for RTBF. Interestingly, the accuracy of the ellipticity is no longer influenced by the absolute level of incoherent noise but by the difference between the incoherent noise ratio on vertical and horizontal components. A method based on the second derivative of the beam power versus the radial wavenumber is finally proposed to experimentally measure the noise ratio difference, which allows experimental values to be corrected. The methodology is compared with classical vertical beamforming and RTBF for a synthetic case and three experimental data sets.

SummaryThe western Yunnan is located in the SE Tibetan Plateau, and is characterized by the active Tengchong volcano (TCV), complex crust-mantle coupling and intense earthquakes. To elucidate tectonism in the western Yunnan, we construct a 3D S-wave velocity model to 80 km depth via ambient noise tomography using dense seismic stations. Our model shows significant low-velocity anomalies at different depths in the crust and uppermost mantle. Compared with the results of previous regional tomography, we image low-velocity anomalies consistent with a large-scale source of partial melts in the uppermost mantle beneath the Tengchong and Baoshan blocks, rather than just below the Tengchong block. Our results also reveal a magma chamber extending from the shallow subsurface to the lower crust beneath the TCV, which is fed by the mantle source. Based on these findings, we propose that the mantle source and crustal magma chamber form a multi-scale magma system. Moreover, the mantle source is potentially resulted from asthenospheric upwelling, which is related to the subduction of the Indian slab. In addition, our model shows that the 1976 M7.4 and M7.3 Longling earthquakes occurred near a magma chamber. Thus, fluids from the magma chamber likely reduced the frictional coefficient on the seismogenic fault and caused the Longling earthquakes.

SummaryDistinguishing between different types of seismic events is a task typically performed manually by expert analysts and can thus be both time- and resource expensive. Analysts at the Swedish National Seismic Network (SNSN) use four different event types in the routine analysis: natural (tectonic) earthquakes, blasts (e.g. from mines, quarries and construction) and two different types of mining-induced events associated with large, underground mines. In order to aid manual event classification and to classify automatic event definitions, we have used fully-connected neural networks to implement classification models which distinguish between the four event types. For each event, we band-pass filter the waveform data in twenty narrow frequency bands before dividing each component into four non-overlapping time windows, corresponding to the P-phase, P-coda, S-phase and S-coda. In each window we compute the root-mean-square amplitude and the resulting array of amplitudes is then used as the neural network inputs. We compare results achieved using a station-specific approach, where individual models are trained for each seismic station, to a regional approach where a single model is trained for the whole study area. An extension of the models, which distinguishes spurious phase associations from real seismic events in automatic event definitions, has also been implemented. When applying our models to evaluation data distinguishing between earthquakes and blasts we achieve an accuracy of about 98% for automatic events and 99% for manually analyzed events. In areas located close to large underground mines, where all four event types are observed, the corresponding accuracy is about 90% and 96%, respectively. The accuracy when distinguishing spurious events from real seismic events is about 95%. We find that the majority of erroneous classifications can be traced back to uncertainties in automatic phase picks and location estimates. The models are already in use at the SNSN, both for preliminary type predictions of automatic events and for reviewing manually analyzed events.

SummaryIn this study, we analyze the 2020 seismic swarm that lasted two months and occurred between the Tancítaro and the Paricutin volcanoes in the Michoacán Guanajuato Volcanic Field, Mexico. We developed a new method to automatically detect and locate about 100,000 earthquakes, enabling us to track the magma migration through narrow dykes. Additionally, we reveal the presence of two magma reservoirs from two seismic noise tomography results. The first reservoir is located from 8 km to 20 km below sea level and beneath the Tancítaro volcano and probably corresponds to a complex network of dykes and sills. This crustal reservoir is fed by a mantle reservoir with a wide horizontal extension between 35 and 50 km below sea level. The seismic swarm initiated beneath the Tancítaro summit in the lower portion of the crustal magma reservoir. At this stage, the seismicity migration was mainly horizontal, which we interpret as its response to the higher normal stress caused by the gravitational load of Tancítaro. Once the magma was displaced laterally from beneath Tancítaro, magma migration became more vertical. The swarm reached the upper portion of the crustal magma reservoir but did not escape it. We also reveal the effect of a distant but strong tectonic earthquake on the seismic swarm. Before its occurrence, magma migration followed several paths; afterwards, it became more focused along a single path. Finally, after the swarm, we observed a second type of seismicity called post-swarm seismicity, with a lower earthquake rate but with higher magnitudes. The hypocenters were diffuse and horizontally centered on the previous swarm location. Furthermore, some earthquakes were aligned along shallow faults, generating a high seismic risk to the different Tancítaro nearby localities.

SummaryThe correlation between the magnitudes of earthquakes is a scientific question usually investigated by statistical seismologists. In the last decades, two opposite answers have been given to the problem of correlation of magnitudes: there is no correlation, i.e. the magnitude of earthquakes in a seismic catalog can be considered as random sampling from the magnitude distribution; or, there is a correlation, i.e. the magnitude of a seismic event can influence the magnitude of the successive event. Here we used the Southern California seismic catalog, properly treated to remove incompleteness, to answer the question, finding no significant correlation between magnitudes.

SummaryThis paper explores the applicability of Ensemble Kalman Inversion (EKI) with level-set parameterization for solving geophysical inverse problems. In particular, we focus on its extension to induced polarization (IP) data with uncertainty quantification. IP data may provide rich information on characteristics of geological materials due to its sensitivity to characteristics of the pore-grain interface. In many IP studies, different geological units are juxtaposed and the goal is to delineate these units and obtain estimates of unit properties with uncertainty bounds. Conventional inversion of IP data does not resolve well sharp interfaces and tends to reduce and smooth resistivity variations, while not readily providing uncertainty estimates. Recently, it has been shown for DC resistivity that EKI is an efficient solver for inverse problems which provides uncertainty quantification, and its combination with level set parameterization can delineate arbitrary interfaces well. In this contribution, we demonstrate the extension of EKI to IP data using a sequential approach, where the mean field obtained from DC resistivity inversion is used as input for a separate phase angle inversion. We illustrate our workflow using a series of synthetic and field examples. Variations with uncertainty bounds in both DC resistivity and phase angles are recovered by EKI, which provides useful information for hydrogeological site characterization. While phase angles are less well-resolved than DC resistivity, partly due to their smaller range and higher percentage data errors, it complements DC resistivity for site characterization. Overall, EKI with level set parameterization provides a practical approach forward for efficient hydrogeophysical imaging under uncertainty.

SummaryThe fundamental behaviors of rock magnetism (piezomagnetic characteristics) at the field scale have not been confirmed because conventional experiments can be performed only in the laboratory. Here, the periodic extraction and injection of gas in the Hutubi ultralarge underground gas storage (UGS) system are used to simulate the stress loading and unloading process at the field scale. We treat 26 gas wells in the UGS system as a multipoint-source Mogi model and calculate models of the piezomagnetic field generated during the operation of the UGS system. These models show that the local magnetic field (LMF) in the southern and central areas of the UGS system showed positive changes. In contrast, the northern area showed negative changes, and the amplitude of the negative changes was smaller than that of the positive changes. Changes in the Curie point depth and gas volume do not significantly alter the spatial characteristics of the LMF.

SummaryStress drop is a proxy of understanding earthquake source process, and it is controversial whether the stress drops of induced earthquakes associated with hydraulic fracturing and injection activities are similar to those of tectonic earthquakes. The measurement of stress drops is usually biased due to the limitations of observation means, or hidden issues in the estimation approaches. Utilizing a local short-period seismic network, we investigate the stress drops of induced earthquakes in Weiyuan Shale Gas Field in Sichuan Province, China from 2019 to 2020. Totally 11844 earthquakes are involved in the analysis, and their stress drops are obtained using an improved approach on the basis of the traditional spectral decomposition method combined with a global optimization algorithm to avoid stacking of spectra that is found leading to source parameter underestimation. We divide the studied area into 3 subareas, and the results show strong stress drop heterogeneity across the entire region. We obtain an average stress drop of 2.29MPa, piecewise stress drop dependence to earthquake magnitude, and complex depth dependence pattern. Our results indicate that stress drops of induced earthquakes are overall consistent with the induced earthquakes in other areas as well as tectonic earthquakes in different environments. Meanwhile, the complexity in the stress drop dependence to depth possibly reflects the variability of stress drops for different earthquake triggering mechanisms.

SummaryThis 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, e.g. 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, characterised 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 characterisation and RN diffusion study has been initiated to investigate various approaches of heterogeneity characterisation 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 characterised 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.