Spatial Heterogeneity of Pore Structure in the Crustal Section of the Samail Ophiolite: Implications for High VP/VS Anomalies in Subducting Oceanic Crust

Seismic surveys along subduction zones have identified anomalously high ratio of P‐ to S‐wave velocity (VP/VS) in the subducting oceanic crust that are possibly due to the presence of pore water. Such interpretations postulate that the pore structure is homogeneous at the scale of the seismic wavelength. Here we present the first statistical evidence of a heterogeneous pore structure in oceanic crust at scales larger than laboratory samples. The spatial correlation of measured bulk density profiles of the crustal section of the Samail ophiolite suggests that the pore structure is heterogeneous at scales smaller than ∼1 m. Wave‐induced fluid flow cannot follow the loading during the seismic wave propagation at this estimated heterogeneity, which implies that fluid‐filled microscopic pores and cracks have a limited impact on the observed high VP/VS anomalies in the subducting oceanic crust. Large‐scale cracks may therefore play an important role in shaping these anomalies.


Introduction
Seismic surveys in subduction zone environments have detected anomalously high ratios of P-to S-wave velocity (V P /V S ) within the subducting oceanic crust (e.g., Audet & Schwartz, 2013;Audet et al., 2009).These anomalies are attributed to the presence of pore fluids, given the influence of pore fluids on the elastic properties of rocks (Peacock et al., 2011).Such seismic anomalies have been associated with seismic activity in subduction zones, as water plays a dominant role in controlling the frictional and rheological properties of rock (Saffer & Tobin, 2011).
Numerous laboratory experiments have been conducted to identify the factors contributing to the observed seismic velocities in oceanic plates, thereby revealing the relationship between the elastic properties of rock and pore fluids (e.g., Akamatsu, Nagase, et al., 2023;Christeson et al., 2016;Hatakeyama & Katayama, 2020;Violay et al., 2010).The elastic wave velocities that are measured in the laboratory are generally determined using ultrasonic transducers with a frequency of 0.1-1 MHz, which means that the measured values reflect the pore structure at the millimeter scale.However, the seismic signals that are utilized in geophysical explorations are typically in the 1-100 Hz frequency band, which corresponds to wavelengths that are larger than tens of meters and reflect the macroscopic structures of the rock mass.It is therefore essential to bridge the spatiotemporal gap between the seismic-scale macroscopic structure of the rock mass and laboratory-scale pore structures to quantitatively interpret the observed field-based seismic velocity structure in terms of pore/crack microstructure.
of the Samail ophiolite is more spatially heterogeneous than previously assumed • The effect of fluid-saturated microcracks on low-frequency seismic velocities is modeled as an unrelaxed condition for this heterogeneity • The high V P /V S anomaly in the subducting oceanic crust can be explained by both microcracks and large-scale cracks

Supporting Information:
Supporting Information may be found in the online version of this article.
A rock mass is generally assumed to be statistically homogeneous at the macroscopic scale to allow for the application of theoretical models, such as effective medium theory (Kachanov, 1994;O'Connell & Budiansky, 1974) and poroelastic theory (Gassmann, 1951), when interpreting the seismic velocity structure.This means that a representative elementary volume (REV) exists, and a physical property is homogeneous and independent of further increases in volume (Bear, 2013;Guéguen et al., 2006).Based on laboratory results using hand-sized specimens, conventional models and interpretations of seismic velocity assume that the REV is on the millimeter scale, even when considering the scale of seismic wavelengths.However, many laboratory studies that have analyzed discrete core samples from the seafloor and ophiolites have commonly detected a wide range of physical properties and porosity, even for samples that were collected from similar depths and lithologies (Akamatsu, Nagase, et al., 2023;Gilbert & Salisbury, 2011;Karato, 1983;Katayama et al., 2021;Violay et al., 2010).This may suggest that the microscopic crack/pore distribution within oceanic crust is actually heterogeneous.It is therefore unclear whether the pore structure at the laboratory-scale REV can be used to model the macroscopic seismic velocity structure.
Here we estimate the scale of the macroscopic REV of the microcrack/pore distribution in oceanic crust using the bulk density profile from core samples collected through the crustal section of the Samail ophiolite during the Oman Drilling Project.The bulk densities of the core samples were measured at a 4-cm interval through ∼400-mlong core using a whole-round multi-sensor core logger (MSCL-W) system (Kelemen et al., 2020b).We quantify the heterogeneity of the microscopic pore structure over the macroscopic scale by using the REV, based on the depth variations in the bulk density that reflect the spatial autocorrelation of pore structure.We then investigate the validity of theoretical models, such as effective medium theory and poroelastic theory, using this estimated REV to provide an interpretation of the imaged low-frequency seismic velocity structure.We propose that the observed high V P /V S anomaly in subducting oceanic crust cannot be explained by microcracks alone, as largescale high-permeability cracks also play an important role.

Samples and Data
The International Continental Drilling Program (ICDP) Oman Drilling Project obtained continuous core samples of the crustal section and upper mantle section of the Samail ophiolite, which is one of the best examples of paleooceanic lithosphere in the world.The recovery rate of the drillcore was ∼100% (Kelemen et al., 2020a).Here we focus on Hole GT3A, which was drilled through the transition zone from sheeted dikes to underlying gabbros in the crustal section of the Samail ophiolite.Basalt, diabase, and gabbro are the main lithologies in the ∼406-m core section (Figure 1).The recovered cores were systematically described and analyzed in the laboratory onboard the drilling vessel Chikyu (Kelemen et al., 2020a).The recovered core samples were carefully reconstructed, shrinkwrapped with core liners, and run through a whole-round multi-sensor core logger (MSCL-W, GEOTEK Ltd.) to measure their physical properties before the core was split and cut into discrete cubic samples.The MSCL-W system automatically measures various physical properties of the core, including the magnetic susceptibility, natural gamma rays, resistivity, P-wave velocity, and bulk density, at an interval of ∼4 cm.We focused on the bulk density in this study, as it is mainly dependent on the porosity and crack density.Note that complete welllogging data that can be compared to our results, such as neutron porosity and sonic velocities, were not collected in this borehole.
The bulk density of the drillcores was determined by gamma ray attenuation (Kelemen et al., 2020a).The gamma ray beam, which is produced by a 137 Cs gamma ray source at a radiation level of 370 MBq radiation level within a lead shield with a 5-mm collimator, is directed through whole-round cores.The gamma ray detector includes a scintillator and an integral photomultiplier tube to record the gamma rays that pass through the whole-round core.The bulk density (ρ b ) is calculated as follows: where I 0 is the gamma ray source intensity, I is the measured gamma ray intensity passing through the sample, μ is the Compton attenuation coefficient, and d is the sample diameter.μ and I 0 are treated as constants, such that ρ b can be calculated from I. The gamma ray detector is calibrated with a sealed calibration core, which is a standard core liner that is filled with pure water and aluminum cylinders of various diameters.The bulk density measurements on the core samples are conducted every 4 cm for 4 s, with a spatial resolution of 5 mm. Figure 1 shows

Geophysical Research Letters
10.1029/2023GL106943 the depth variations in the measured bulk density of Hole GT3A.The measured values generally range from 2.3 to 3.2 g/cm 3 , with some unusually low values that are less than ∼2 g/cm 3 .These values can be attributed to low gamma ray attenuation due to missing core sections or fragmentation of the cores, and do not reflect the true density.These low values are therefore considered outliers.
The remaining data may also contain measurement errors and/or noise, which can obscure the true spatial variations in bulk density.We used Markov Random Field (MRF) modeling, which is a kind of Bayesian estimation, to overcome these uncertainties by reconstructing the bulk density profile.The MRF model allows for data-driven noise removal and accurate reconstruction of the true data by assuming continuity in the spatial distribution of the physical quantity (Geman & Geman, 1984), and has been widely applied to various fields in the natural sciences, including geosciences (Kuwatani et al., 2014).
Figure 1 also shows the depth variation in the bulk density that was estimated via the MRF model.See the Supporting Information S1 for details of the MRF modeling.The reconstructed bulk densities are much smoother than the raw data and should be more accurate, which can be supported by the successful hyperparameter estimation (Figure S1 in Supporting Information S1).The bulk density profiles that were reconstructed via MRF modeling exhibit no systematic changes with depth, although the values tend to be slightly lower in the deepest section of the hole, which may coincide with the lithological transition from sheeted dikes to gabbros.

Estimation of the REV
We evaluated the spatial heterogeneity and REV size of microscopic pore structure using the bulk density profile of Hole GT3A that was reconstructed via MRF modeling.The REV is defined as the minimum volume for which a physical property is homogeneous and independent of further increases in volume (Bear, 2013).Here, we estimate the size of the REV based on the spatial autocorrelation of the bulk density which is obtained via empirical semi-variograms (AlRatrout et al., 2018;Jackson et al., 2020).A semi-variogram is a statistical tool that quantifies how the variance of a variable changes with increasing distance between data points.It is used in geostatistical spatial analysis to measure the degree of dissimilarity or variability in geological characteristics, including the fracture density of a rock and a borehole (Koike & Ichikawa, 2006).The semi-variogram for one-dimensional data is calculated as follows: where N(h) is the number of measurement point pairs that are separated by a distance h (i.e., the lag distance), and x i and x i+h are the physical quantities at measurement points i and i + h, respectively.Here we calculate the semivariograms through the entire core with the maximum lag distance set to 100 m and with 20 m of overlap between the calculated intervals.
The semi-variograms calculated from the bulk density profile that was reconstructed via the MRF model are shown in Figure 2. The semi-variograms initially increases with increasing lag distance, which indicates that the bulk densities are spatially correlated within those length scales.The semi-variograms then approach plateaus, which represent the variance at the largest scale that is considered.These plateaus indicate the disappearance of the spatial correlation of the bulk density, and suggest statistical homogeneity.The distances to the plateaus can be considered a lower bound of the REV length (Jackson et al., 2020).The distances to the plateaus of the bulk density of Hole GT3A are 0.5-1.4m, which suggests that the effect of microscopic pores/cracks on the physical properties, at the scale of ∼100 m, can be modeled by the pore structure at the meter scale.Nevertheless, it should

Geophysical Research Letters
10.1029/2023GL106943 be noted that the macroscopic lithological variations observed in Hole GT3A (Figure 1) may be indicative of heterogeneity at a scale larger than this borehole.

Wave-Induced Fluid Flow (WIFF) Within the REV
The size of the REV can matter when modeling the elastic properties of a rock containing fluid-saturated pores/cracks.When porous/cracked rocks are fluid-saturated, the loading associated with the elastic wave propagation induces a pore-pressure gradient in the pore network (Mavko et al., 2020).The pore pressure is diffused by fluid flow, which is referred to as WIFF.The effective elasticity of fluid-saturated rocks depends on the rate of loading of the REV and the rate of pore-pressure diffusion across the REV (Guéguen & Kachanov, 2011).This means that the elastic wave velocity depends on both the frequency and the spatial scale of the REV.If the fluid flow has no time to take place within a given REV during wave propagation, then the pore pressure within the REV is different among the cracks.The measured elastic properties of this rock are referred to as "unrelaxed" moduli in this case.These moduli can be modeled via effective medium theory, which postulates non-uniformity of the pore pressure within the REV (e.g., Kachanov, 1994).However, if there is enough time for fluid flow to occur between cracks during wave propagation, then the pore pressure can be equilibrated on the scale of the REV.
The measured elastic properties of the rock are referred to as "relaxed" moduli in this case, and are within the validity range of the concept of poroelastic theory.Nevertheless, there is no time for fluid flow to extend beyond the scale of the REV, and the pore pressure between REVs varies at the macroscopic scale.This situation corresponds to the "undrained condition" of poroelastic theory, and the undrained elastic moduli can be predicted by Gassmann's equations (Gassmann, 1951).The undrained moduli are typically smaller than the unrelaxed moduli, which means that the effect of water-saturated pores/ cracks on the seismic velocity is enhanced under undrained conditions.
We calculated the time required for fluid flow to diffuse over the REV length that was estimated for the microscopic pore structure of Hole GT3A to choose an appropriate model for predicting the corresponding seismic velocity in oceanic crust.We then considered whether the pore pressure can be equilibrated within the REV (i.e., undrained or unrelaxed) for typical frequency ranges.The time required for the fluid to diffuse over distance L within a rock containing cracks and pores is calculated as follows: where η is the fluid viscosity, K d is the drained bulk modulus, and k is the permeability (Guéguen & Palciauskas, 1994).The critical frequency at which the undrained/unrelaxed transition f c occurs is the inverse of Equation 3: If the frequency of the elastic wave is greater than f c , then the fluid diffusion over length L is slower than the rate of loading, resulting in unrelaxed conditions.Conversely, if the frequency is less than f c , then the fluid diffusion is faster than the rate of loading, resulting in relaxed but undrained conditions.
We used the measured elastic and electrical properties of discrete cubic core samples from Hole GT3A (Akamatsu, Nagase, et al., 2023) to obtain the k and K d values for Hole GT3A.Crack percolation theory states that the permeability of rock containing randomly oriented cracks can be defined as follows: where w is the crack aperture and F is the electrical formation factor, which is a parameter that represents the geometry of the crack network (Guéguen & Dienes, 1989).We used the F values derived from the electrical resistivity of cubic samples under brine-saturated conditions.We assumed a 1-μm crack aperture, which yielded k = 5.1 × 10 17 m 2 to 6.5 × 10 16 m 2 for Hole GT3A, with a median value of 1.4 × 10 16 m 2 .K d was set to 13.6 GPa, which was derived from the P-and S-wave velocities under dry conditions.η was assumed to be 10 4 Pa.s, which is plausible for in situ crustal conditions (Pimienta et al., 2018).
The critical frequencies, which were calculated from Equation 4and the above-mentioned values, are shown in Figure 3 (dark shaded region) as a function of length.The ranges of typical micro-scale REVs (10 3 to 10 2 m) and the macro-scale REVs that were determined in this study (0.5-1.4 m) are also illustrated, along with the typical frequency ranges of laboratory ultrasonic (10 5 -10 6 Hz), borehole sonic logging (10 3 -10 4 Hz), and seismic surveys (10 0 -10 2 Hz).The pore pressure within the REV is unrelaxed in the case of the conditions in the region above the critical frequencies where effective medium theory can be applied, whereas it is relaxed (but undrained) in the case if the conditions in the region below the critical frequencies where poroelastic theory can be applied.The ultrasonic velocity reflects unrelaxed conditions at the conventional microscopic REV scale, whereas the seismic velocity reflects undrained conditions (dashed blue rectangle in Figure 3).The sonic logging velocity coincides with the transition between unrelaxed and undrained conditions.However, our analyses suggest that the REV of the microscopic pore structure at the macroscopic scale can be up to ∼1 m in the oceanic crust (red region in Figure 3).This estimate indicates that the pore pressure conditions within the REV for microscopic pores structure in oceanic crust for seismic, sonic, and laboratory ultrasonic velocities might all be unrelaxed.Achieving relaxed conditions for seismic frequency would require a permeability of at least ∼10 14 m 2 .Therefore, effective medium theory can be used to model the effect of microscopic porosity (i.e., matrix porosity) on seismic velocity.

Implications for High V P /V S Anomalies in the Subducting Oceanic Crust
Seismic surveys have revealed anomalously high V P /V S values of ∼2.4 within the low-velocity zones in the interior of subducting oceanic crust (Audet & Schwartz, 2013;Audet et al., 2009), even though the pressures at the depths of these low-velocity zones are high enough to result in crack closure.These anomalies possibly indicate that the pore pressure is close to the lithostatic pressure, which would be high enough to keep the cracks open.Since the pore pressure largely controls the rheological properties of rock, it has been suggested that such high V P /V S values are linked to seismic activity along subduction zones (Saffer & Tobin, 2011).The observed high V P /V S anomalies have been related to highly microfractured rock under undrained conditions, based on laboratory measurements of frequency-dependent elastic moduli (Pimienta et al., 2018).However, our analysis implies that the seismic velocity of oceanic crust should be modeled under unrelaxed conditions, given the estimated heterogeneity of pore structure of Hole GT3A.
Figure 4 shows the V P /V S as function of crack density under unrelaxed and undrained conditions, as calculated using differential effective medium theory and poroelastic theory, respectively (Le Ravalec & Guéguen, 1996).
Here the crack density represents the degree of fracturing in a given REV.The mean crack aspect ratios are set to typical values (10 3 and 10 2 ) in the calculations, and the crack-free P-and S-wave velocities are assumed to be 7.0 km/s and 3.8 km/s, respectively (Akamatsu, Nagase, et al., 2023).V P /V S increases rapidly with crack density for both aspect ratios under undrained conditions, thereby explaining the high V P /V S values of subducting oceanic crust (Audet & Schwartz, 2013).This result is in agreement with the interpretation based on the laboratory measurements of low-frequency elastic properties under the assumption of a small REV (Pimienta et al., 2018).
Figure 3. Critical frequencies at which the unrelaxed/undrained transition occurs in rock as a function of length based on the observed physical properties of Hole GT3A (see the main text for details).The length ranges of the macroscopic representative elementary volume (REV) that are determined in this study and the conventional microscopic REV are highlighted.The typical frequency ranges for laboratory ultrasonic measurements (10 5 -10 6 Hz), borehole sonic logging (10 3 -10 4 Hz), and seismic surveys (10 0 -10 2 Hz) are also shown along the right vertical axis.

10.1029/2023GL106943
Conversely, V P /V S increases relatively slowly with crack density under unrelaxed conditions, such that either an extremely high crack density or a small aspect ratio are required to explain the observed anomalies.
Here we compare these model predictions with data from discrete core samples from Hole GT3A, which reflect unrelaxed conditions (Akamatsu, Nagase, et al., 2023).The measured values on the discrete core samples do not reach the observed V P /V S anomalies in the subducting crust.This discrepancy suggests that the observed anomalies cannot be explained solely by the identified microcrack characteristics of Hole GT3A.Otherwise, the structure of the microcrack network associated with brittle fracturing in subducting crust could be largely homogeneous and characterized by a microscopic REV.
The other explanation for this gap between the seismologically observed and laboratory measured V P /V S values is the occurrence of large-scale cracks.
Here we focused on the heterogeneity of the microscopic pore structure within Hole GT3A.However, natural geological systems likely contain multiple scales of porosity structures, as opposed to just microscopic pores and cracks (Bailly et al., 2019;Ringrose et al., 2008).The core samples from Hole GT3A contain numerous millimeter-wide mineral veins, which could be interpreted as cracks that once acted as fluid pathways within the oceanic crust (Akamatsu, Katayama, et al., 2023).Such cracks may have much higher permeabilities of up to ∼10 10 m 2 (Akamatsu, Katayama, et al., 2023;Hosono et al., 2022;Ioannou & Spooner, 2007), as permeability is strongly dependent on the crack aperture (Guéguen & Palciauskas, 1994).WIFF can rapidly spread within a given REV that possesses such a high permeability, even though both the scale of heterogeneity and REV of the large-scale cracks might be a few orders of magnitude larger than those of microcracks (Bailly et al., 2019;Koike & Ichikawa, 2006).While the pore-pressure condition for microscopic crack/pore networks remains unrelaxed, that for the large-scale crack networks could be undrained for relatively low-frequency sonic/seismic velocities if the REV length of the large-scale cracks is less than several tens of meters.This suggests that the effect of large-scale cracks can be greater than that of microcracks, thereby leading to the anomalously high V P /V S values in subducting oceanic crust.Such large-scale cracks may play an important role for subduction zone seismicity that are associated with high pore pressure and fluid migration (Kodaira et al., 2004;Saffer & Tobin, 2011;Shelly et al., 2006).

Conclusions
Spatial variations in the bulk density of the drillcore from the crustal section of the Samail ophiolite, Oman, were used to estimate the REV of the microcrack distribution in oceanic crust.The estimated REV length is 0.5-1.4m, which is much larger than that assumed by conventional effective medium modeling of seismic velocity.The effect of WIFF on low-frequency elastic moduli is characterized as an "unrelaxed condition" at the estimated REV scale.This implies that fluid-filled microcracks have a limited impact on the observed high V P /V S anomalies in subducting oceanic crust and that large-scale cracks potentially have a larger impact.Exploring the scale dependence of seismic velocity may be useful for understanding the temporal seismic activity that is associated with fluid migration through cracks in subduction zones.(Kachanov, 1994) and poroelastic theory (Gassmann, 1951), respectively.The numbers in parentheses are the mean crack aspect ratios.The gray area represents the observed values in subducting oceanic crust (Audet & Schwartz, 2013).The circles denote measured values for discrete core samples from Hole GT3A, which reflect unrelaxed conditions (Akamatsu, Nagase, et al., 2023).

Figure 1 .
Figure 1.Depth variations in the bulk density and borehole stratigraphy of Hole GT3A.The gray circles are the raw data, which were measured with the MSCL-W system, and the blue lines are the values that have been denoised using the Markov Random Field model.

Figure 2 .
Figure 2. Semi-variograms of the bulk density profile of Hole GT3A, reconstructed the Markov Random Field model.Darker colors represent greater depths.The gray region indicates the distances to the plateaus in the semi-variograms, which can be considered a lower bound of the representative elementary volume length. k

Figure 4 .
Figure 4. V P /V S as a function of the crack density parameter under unrelaxed and undrained conditions, modeled via differential effective medium theory(Kachanov, 1994) and poroelastic theory(Gassmann, 1951), respectively.The numbers in parentheses are the mean crack aspect ratios.The gray area represents the observed values in subducting oceanic crust(Audet & Schwartz, 2013).The circles denote measured values for discrete core samples from Hole GT3A, which reflect unrelaxed conditions(Akamatsu, Nagase, et al., 2023).