Fracture Shape and Orientation Contributions to P-Wave Velocity and Anisotropy of Alpine Fault Mylonites

1. Introduction

The Alpine Fault is located along the West Coast of the South Island, New Zealand. It marks the transpressional plate boundary between the Australian and the Pacific Plates (Figure 1). The fault hosts many earthquakes a year, but a large earthquake (M_w > 7.0) occurs every <300 years on average with the last one occurring in 1717 AD (Sutherland et al., 2007; Howarth et al., 2018). The majority of the rocks on the Pacific Plate (hanging wall) are schist and fault rocks. Rocks sheared during the collision of the two plates form a series of fault rocks close to the principal slip zone consisting of fault gouge, cataclasite, and variable grades of mylonites which vary according to fault-perpendicular distance. Mylonites do not outcrop in the footwall, but the hanging wall mylonite zone is approximately 1 km thick (Norris and Cooper, 1997). The mylonites are derived from the Alpine Schist and developed foliation parallel and sub-parallel to the fault plane (Toy et al., 2008). In general, the mylonites and schist are composed of quartz, plagioclase, biotite, muscovite, and accessory minerals such as chlorite, garnet, and calcite (Grapes and Watanabe, 1994; Little et al., 2002; Toy et al., 2015; Boulton et al., 2017).

Figure 1. Simplified geological map of the study area and cross-section showing typical shallow-depth sequence of the Alpine Fault rocks (Norris and Cooper, 2007; Schuck et al., 2020). Alpine Fault trace shown in orange; Stony Creek and Smithy Creek are indicated by black triangles. Map is based on the geological map of New Zealand's South Island (Edbrooke et al., 2015) and the digital elevation model of Columbus et al. (2011).

Field passive and active source seismic data has been used to investigate the Alpine Fault geometry and rock physical properties. Such data show a low-velocity zone in the hanging wall that lies parallel to the Alpine Fault (Smith et al., 1995; Stern et al., 2001, 2007; Feenstra et al., 2016; Lay et al., 2016). Within the first 10 km depth, the P-wave velocity of the low-velocity zone is 5.5–6.0 km/s compared to the velocity outside the zone of 6.0–6.5 km/s. This low-velocity zone is commonly interpreted to be due to the presence of fluid-filled microfractures, but experimental and model data needed to validate such a claim are rare (Allen et al., 2017; Adam et al., 2020). The fault damage zone has been investigated using fault guided (trapped) seismic wave observation and modeling (Eccles et al., 2015; Kelly et al., 2017). However, these field studies that aim to characterize the fault geometry and physical properties do not account for seismic anisotropy. It is well-known that assuming isotropic seismic velocity models for anisotropic rocks will result in erroneous modeling and interpretations (Chapman and Pratt, 1992; Godfrey et al., 2002; Tsvankin et al., 2010). It is therefore critical to quantify how the physical properties of Alpine Fault rocks influence elastic wave anisotropy.

Previous studies using seismic data (Savage et al., 2007; Karalliyadda and Savage, 2013), core measurements (Okaya et al., 1995; Godfrey et al., 2000; Christensen and Okaya, 2007; Allen et al., 2017; Simpson et al., 2019; Adam et al., 2020; Jeppson and Tobin, 2020) and numerical modeling (Godfrey et al., 2002; Dempsey et al., 2011; Adam et al., 2020; Simpson et al., 2020) demonstrate that the rock in the vicinity of the Alpine Fault is strongly anisotropic. The anisotropy is attributed to many factors such as mineral crystallographic preferred orientation (CPO), grain shape preferred orientation (SPO), foliation and fractures. The Alpine Fault rocks contain significant (>15 vol%) mica with preferred orientation (Dempsey et al., 2011; Adam et al., 2020). Dempsey et al. (2011) has shown how the volume of aligned mica influences P-wave anisotropy in mylonite from the Alpine Fault. Moreover, grain size insensitive creep-accommodated ductile shear resulted in a CPO of quartz and plagioclase as well as their SPO (Toy et al., 2008). Velocity measurements under confining pressure also show directionally, non-uniform changes in P-wave velocities with depth which implies that aligned fractures are present (Adam et al., 2020; Simpson et al., 2020).

Fractures herein are described as flat-shaped pores. Fractures can form within the rock as intergranular spaces between flat-shaped grains in sedimentary rocks (Loucks et al., 2009), or can form as a secondary structures when the rocks are exposed to differential stress (e.g., Pischiutta et al., 2015). In many tectonic scenarios, fractures can develop as an aligned set. These aligned fractures affect the seismic anisotropy of rocks. The fractures can induce seismic anisotropy in an isotropic host rock, or further enhance pre-existing anisotropy (Kern et al., 2008; Novitsky et al., 2018). The degree of fracture-induced seismic anisotropy enhancement depends on many factors such as fracture volume (porosity), orientation, and shape (Hong-Bing et al., 2013; Vasin et al., 2013). These factors should be incorporated into seismic velocity modeling.

Adam et al. (2020) showed that open grain boundaries and microfractures can enhance the seismic wave anisotropy at all depths of the seismogenic zone in the Alpine Fault. Simpson et al. (2020) predict the P-wave anisotropy of shallow Alpine Fault rocks. However, those studies did not systematically consider the effect of mica CPO and shape/orientation of fractures on seismic anisotropy.

Here we quantify fracture shape, mineral phase, and crystal orientations from synchrotron X-ray microtomography data and electron back-scattered diffraction (EBSD) data. Based on these measurements, we model P-wave speeds and anisotropy for a range of Alpine Fault rocks. The effect of fractures on Alpine Fault rocks is modeled using two numerical methods: a wave propagation and a differential effective medium modeling approach. The methods differ in being dynamic and static respectively, as well as in the way that fractures are in 2D and 3D. Our study compares the models and their limitations for the inclusion of mica CPO and distribution, and fractures to predict seismic wave anisotropy.

2. Materials and Methods

2.1. Electron Back-Scattered Diffraction (EBSD)

Rock samples were collected from Central Alpine Fault outcrops which span shear zone mylonites and schist. The schist sample was collected at Smithy Creek and is part of the Alpine Schist tectonostratigraphic unit (Cooper and Palin, 2018). The protomylonite, mylonite and ultramylonite samples are from Stony Creek (Figure 1). Thin sections of 35 μm thickness are prepared for all samples and subjected to fine SYTON polish to be examined using EBSD. Thin sections for the samples are created for two angles to foliation: 90^o (perpendicular to foliation) and 30^o. No pores or fractures are visible in the specimens under optical microscope or in microphotographs.

EBSD experiments are performed to produce a phase map containing mineralogical and crystallographic orientation data. EBSD data were acquired on 2 by 2 mm areas on the thin sections using a Zeiss Sigma VP FE-SEM fitted with an Oxford Instruments HKL INCA Premium Synergy Inte30 grated ED/BSD system, located at the Otago Microscopy and Nano-Imaging (OMNI). The SEM is operated with an accelerating voltage of 30 keV, an aperture size of 300 μm, and a working distance of approximately 30 mm. The EBSD data were acquired at a 2 μm step size then used to quantify mineral phase composition and crystallographic preferred orientation (CPO).

The EBSD data were processed to reduce noise and fill any data gaps using the MTEX toolbox (Mainprice et al., 2015). First, grain boundaries were reconstructed by grouping same-phase mineral data with less than 10^o misorientation angle. Then, all isolated 1-pixel grains were deleted to eliminate wild pixels (i.e., noise). Grain boundaries were reconstructed again using the same criteria as the previous step. Grains with less than 30% indexed fraction were removed. The indexed fraction is a ratio between the area of the indexed pixel in the grain to the whole grain area. Grains smaller than 10 pixels were also removed. Lastly, in-grain pixels were smoothed and filled based on their surrounding pixels using a spline filter.

Fine-grained phyllosilicate minerals, dominant in mylonite rocks, are mostly unable to be indexed by EBSD and appear in the maps as non-indexed phase (Prior and Mariani, 2009; Dempsey et al., 2011; Adam et al., 2020). To improve indexing of mica, we performed EBSD on an additional cut to the samples at 30^o from foliation with the hope to gain surface area on the mica to improve indexing. Unfortunately, we had about the same volume of non-indexed phase as for the 90^o cut to foliation. Toy et al. (2008) and Dempsey et al. (2011) show the smallest phases in the Alpine Fault mylonites are dominantly mica. We therefore assign the non-indexed phases to mica. Moreover, based on photomicrographs and synchrotron micro-CT imaging, mica within Alpine Fault mylonite (Toy et al., 2008; Adam et al., 2020) is parallel and subparallel to foliation. We therefore, assume that most non-indexed phases (i.e., mica) are parallel to foliation and account for the sub-parallel mica in terms of randomly oriented mica. Microfractures and open grain boundaries would also show as non-indexed phases, however, these are mostly not resolved in the EBSD imaging. It is important to understand that the EBSD data does not have information regarding pores or fractures of these rocks. Since all non-indexed pixels are substituted with mica in our models, the effect of microfractures is later included separately in the finite difference and the wave propagation models.

It is clear that not all micas are perfectly parallel to foliation. To determine the mica orientation, we use the P-wave high pressure data from Adam et al. (2020) where most of the fractures are assumed to be closed with a 0.5% modeled porosity to estimate the volumetric contribution of parallel and randomly oriented mica basal plane orientation. We compare MTEX models of fast and slow P-wave velocities based on the combinations of perfectly and randomly aligned fractures volume to the laboratory data (Supplementary Figure 4). We show that model and experimental data match for a 80:20 volume of parallel:random fractures. We therefore, use this volumetric distribution of mica for all of our models.

2.2. Synchrotron X-Ray Micro-Tomography (Micro-CT)

Micro-CT is a non-destructive method used to investigate the 3D internal structure of an object. The shape of rock microfractures can be extracted and these can statistically be analyzed and converted into mathematical shapes. We use such information to model fracture effects on seismic anisotropy with the wave propagation and differential effective medium models.

Micro-CT data used in this study was collected using beamline 20XU at SPring-8, Japan on a 3 mm height and 1 mm diameter cylinder. The sample was mounted on a rotary stage. X-ray beam of 20 keV energy was shone through the sample onto the scintillator which converts X-ray to visible light. The image was then captured by the CMOS camera detector. The process was repeated over 180^o rotation with 0.1^o step size. The sample cross-sections were reconstructed using a convolution back-projection method, producing a stack of gray-scale which then re-scaled to 8-bit gray-scale images. The data resolution is 0.524 μm.

The micro-CT data were processed using the Avizo software. The gray-scale micro-CT images were cropped into an approximately 250-micron cube. The data were smoothed with a first-degree median filter to reduce noise and artifacts. The filter replaces the CT number in a voxel with the median of the 27 neighboring voxels (voxels within a 1-voxel radius of the center voxel). All pores which are smaller than 8 voxels were also removed, and a hole-filling algorithm was used. After that, pore segmentation was performed using a threshold method based on the CT number. Lastly, 3D pores were reconstructed, and the pore dimension and volume fraction were measured. EBSD and micro-CT data were used to model wave speeds in pore-free and fractured mylonites and schist. These data were combined to compare two modeling approaches: wave propagation and GassDEM modeling.

2.3. Wave Propagation Modeling (EWAVE)

Wave propagation modeling is performed using EWAVE Matlab code (Zhong and Frehner, 2018). The code simulates wave propagation through the EBSD data using an explicit finite-difference time-marching model. The model discretizes the input EBSD data into cells and an elastic stiffness tensor (Table 1) is assigned to each cell based on mineralogy and crystallographic orientation. A Ricker wavelet is digitally propagated across the sample by solving a momentum balance equation and Hooke's law. The wave travel time is recorded and the wave velocity is calculated from the width of the EBSD data and the travel time. The fast and slow velocities are determined by propagating waves parallel and perpendicular to rock foliation, respectively. A detailed explanation is provided in Zhong and Frehner (2018). A time snapshot of wave propagation through a medium is shown in Figure 2. The buffer length and wavelength used in this study are 300 and 100 μm, respectively. P-wave anisotropy is calculated as a ratio between the difference between fast and slow velocity, and the average of the two velocities.

Table 1. Elasticity (in GPa) and density (in kg/m³) of mineral phases present in our samples.

Figure 2. EWAVE model component: EBSD data (“Rock”) is placed between two buffers. A Ricker wave (pink and blue) is propagated from left to right of the model. Two receivers (left: black line, right: red line) are used to record the wave signal and calculate wave travel times. Displacement is dimensionless and does not affect the velocity.

Two types of modeling are performed: (1) a porosity-free EBSD model and (2) the same model but with added fractures defined by specific aspect ratios. EWAVE models EBSD data and could include imaged natural fractures (non-indexed phases) as it propagates a wave through the model (Zhong and Frehner, 2018). However, for our mylonitic rocks, natural fractures were not imaged via EBSD so are not immediately accounted for. We include pore information by extracting natural 3D fracture shapes from micro-CT data. This fracture set is statistically and mathematically simplified into 2D fractures which are deterministically included into the model. These pores are added randomly to the EBSD model at specific aspect ratios parallel to foliation.

2.4. Differential Effective Medium Modeling (MTEX and GassDEM)

A different approach to model the pore-free rock elasticity is by applying the Hill effective medium model implemented in the Matlab MTEX toolbox (Mainprice et al., 2015). The calculation uses the mineral volume fraction and CPO from the EBSD data, together with the single-crystal stiffness tensor of each mineral (Table 1) to predict the whole rock effective stiffness tensor from which wave speeds with direction can be estimated.

To model the effect of pores and fluids on the wave velocities, the pore-free effective stiffness tensor output from MTEX is set as a background rock for GassDEM (Gassmann and DEM) modeling (Kim et al., 2019; Simpson et al., 2020). GassDEM modeling adds fixed aspect ratio inclusions (3D fractures in this case) into the background medium. The process is iterative with additions of small volumes of pores until the target sample porosity is reached (Bruggeman, 1935). Gassmann's equation allows the effect of fluids on the elasticity of the rock to be modeled (Gassmann, 1951; Smith et al., 2003).

We study the effect of fractures on P-wave velocity and anisotropy in two ways: (1) fractures are aligned to foliation and porosity volume is varied and (2) the porosity volume is constant, but the contribution of the fracture orientation is changed for different combinations of aligned and randomly oriented fractures.

3. Results

3.1. Electron Back-Scattered Diffraction (EBSD)

EBSD data show that the shear-zone Alpine Fault rocks are composed of quartz, plagioclase, garnet and a considerable amount of non-indexed phase (Figure 3, Table 2), assumed to be biotite mica as the phyllosilicate minerals in the Alpine Fault mylonites are predominantly biotite (Toy et al., 2015). There is no significant mineralogical difference between the 90^o and 30^o sections nor an improvement of imaging the micas. 80% of the non-indexed mica is assumed to align with foliation and 20% of the non-indexed mica is assumed to be randomly oriented (see section 1.3 in Supplementary Material).

Figure 3. Processed EBSD data of mineral phases distribution of (A,B) schist, (C,D) protomylonite, (E,F) mylonite, and (G,H) ultramylonite cut 90^o and 30^o to foliation, respectively.

Table 2. Mineral composition (vol.%) and mineral density (kg/m³).

3.2. Synchrotron X-Ray Micro-Tomography (Micro-CT)

Micro-CT data of ultramylonite, mylonite, and protomylonite samples are analyzed to quantify porosity and pore aspect ratio. After pore segmentation is performed (Figure 4), the number of pore voxels is divided by the total number of voxels to calculate porosity. The calculated porosity of ultramylonite, mylonite, and protomylonite are 3.36, 2.08, and 4.88%, respectively. A 5% porosity will be used in the following velocity modeling as it covers all the calculated porosity from the micro-CT data. Pore aspect ratio is measured as a ratio of three orthogonal pore axis lengths (e.g., X:Y:Z where X≥Y≥Z). The three numbers represent the longest axis, the second-longest and the shortest axial length, respectively (Figure 5). Two axial ratios are extracted from the segmented pores: elongation and flatness (Figure 6). Elongation is a ratio between the second-longest and the longest axial length (i.e., Y/X) and flatness is a ratio between the shortest to the second-longest axial length (i.e., Z/Y).

Figure 4. (A) An example of pore segmentation (blue) on a 2D micro-CT slice for a mylonite. (B) The connected 3D pore network (pink) in the sample. The size of the cubic sub-volume of the micro-CT data analyzed and shown here is 470 X 400 X 400 μm³.

Figure 5. Ellipsoid examples representing fractures of (A) 5:2:1 and (B) 5:5:1 aspect ratio.

Figure 6. Histograms of pore elongation and flatness in (A,B) ultramylonite, (C,D) mylonite, and (E,F) protomylonite.

The elongation and flatness histogram plots of all three rocks are similar. Elongation plots are right-skewed with mode around 0.16 with the mean and median range from 0.36 to 0.40 and 0.32 to 0.38, respectively. The flatness for the ultramylonite and the protomylonite are almost symmetric with 0.48 mean and median. For the mylonite however, the flatness plot is slightly left-skewed. The mean and median of the mylonite flatness are 0.54 and 0.56, respectively. Pores with small elongation and flatness would represent low aspect ratio fractures, which are thin, long and narrow. Pores with elongation and flatness close to 1 represent close-to-spherical pores. Based on the histogram analysis, the average pore aspect ratio determined from the mean flatness and elongation is 5:2:1. This fracture shape is defined as a flat-shaped pore (fracture) and used for the following P-wave velocity modeling.

3.3. Wave Propagation Modeling (EWAVE)

Fast and slow P-wave velocities are modeled using the EWAVE code. Fracture porosity is added randomly as a 2D projection of the extracted micro-CT fractures of aspect ratio 5:2:1 (Figure 7). Fractures are added with their long-axis parallel to foliation as rectangles of aspect ratio 5:1. For this study, fractures are added by randomly replacing pixels in the EBSD data with air/water. The influence of fracture porosity volume is studied by adding fractures in 1% increments up to a total porosity of 5%. In this modeling, pores are dry (filled with air). For such low porosity, modeling air- or water-filled fractures results in a velocity difference of less than 2% of the air-filled fracture model. The pore-free fast velocities are similar for all rocks at 6.5 km/s, while the slow velocities range from 5.5 to 6.0 km/s depending on the lithology (Figure 8A). As porosity is added to the model, both fast and slow velocities decrease linearly. The fast velocities decrease slightly to 6.3 km/s at 5% porosity. The slow velocities however, decrease rapidly to 4.5 km/s. The P-wave anisotropy is also calculated (Figure 8C), ranging from 10 to 15% and increasing linearly to 30–35% for a 5 % fracture porosity.

Figure 7. (A) Schist EBSD data cut perpendicular to foliation and (B) the EWAVE model of the schist data with 5% fractures included. The color represents the phase density. The air-filled fractures are shown in blue.

Figure 8. Modeled P-wave velocity and anisotropy with (A,C) EWAVE model (5:1 fracture aspect ratio) and (B,D) GassDEM model (5:2:1 fracture aspect ratio). U, M, P, S refer to ultramylonite, mylonite, protomylonite and schist, respectively. -fast and -slow refer to fast and slow velocity.

3.4. Differential Effective Medium Modeling (MTEX and GassDEM)

The fast and slow velocities estimated with MTEX have similar values to the EWAVE modeling at approximately 6.5 and 5.5–6.0 km/s, respectively. Foliation-parallel ellipsoidal fractures of aspect ratio 5:2:1 are added to the model for porosity volume from 0.1 to 5% at 0.1% increments. The fast velocities decrease from 6.5 to 6.3 km/s at 5% porosity (Figure 8B), while the slow velocities decrease from 5.5 to 5 km/s. The pore-free P-wave anisotropy ranges from 10 to 15% and increases to 15–20% for a 5% porosity. Although EWAVE and GassDEM/MTEX use the same fracture volume and orientation, slow velocities - and thus P-wave anisotropy—greatly differ between the models. The P-wave velocity and anisotropy results for the section cut at 30^o from foliation and water-filled samples are presented in the Supplementary Material.

The previous models assumed that all fractures are the same size and align with foliation. Fractures in real rocks however, vary in size and orient in different directions. For this, we use the modified GassDEM code (Simpson et al., 2020) to model how different combinations of foliation, parallel and randomly oriented fractures, as well as fracture size, affect P-wave velocities. Two fracture aspect ratios are used: a 5:2:1 aspect ratio which used in this study, and a 50:50:1 ratio as a representative of flat fractures. Eleven fracture orientation combinations are modeled from 0% aligned (random) fractures, to 100% aligned fractures, with combinations of these. Because all rocks show similar velocity and anisotropy trends, only the protomylonite data are used in this model. The P-wave velocities and anisotropy of all fracture models are calculated as a function of propagation angle to foliation. At 90^o the P-wave propagates parallel to foliation, and at 180^o perpendicular to it.

For the flat-shaped pore of 5:2:1 aspect ratio (Figure 9A), the fracture alignment can increase and decrease the fast and slow P-wave velocity, respectively, by approximately 0.4 km/s. Velocity changes become more prominent for the flatter aspect ratio fractures (50:50:1) resulting in overall slower velocities for all azimuthal directions compared to fractures of 5:2:1 aspect ratio (Figure 9B). The aligned fractures with a 50:50:1 aspect ratio increase the fast velocity by 3 km/s and decrease the slow velocity by 1.5 km/s. Fracture porosity has more influence on the P-wave anisotropy as more fractures are aligned (Figures 9C,D). However, when fracture alignment drops below 30% the anisotropy decreases with porosity instead.

Figure 9. P-wave velocity and anisotropy with angles of propagation to foliation for a 5% porosity protomylonite from GassDEM modeling. Two aspect ratios are compared: (A,C) 5:2:1 and (B,D) 50:50:1. The color represents fracture alignment to foliation, from all fractures being randomly oriented (0%) to perfectly aligned fractures to foliation (100%), and combinations of these. The black dashed lines are results for a rock with 100% spherical pores (1:1:1 aspect ratio) for reference.

The modeled velocity using spherical pores (black dashed lines) is presented for comparison (Figure 9). The spherical pore model yields a similar fast velocity to the model including 30% aligned fractures of 5:2:1 aspect ratio. The slow velocity of the 100% randomly oriented fracture model coincides only with the model that only contains spherical pores. For the 50:50:1 models, the fast velocity from the spherical pore model is comparable to the perfectly aligned 50:50:1 pore aspect ratio model. However, the 50:50:1 pore model's slow velocity is significantly lower than that of the spherical pore model. The anisotropy of the spherical pore model remains almost constant and equal to the pore-free anisotropy because spherical pores reduce P-wave velocity equally, regardless of direction as opposed to fractures. If only 10–20% of the fractures align with foliation, then the P-wave anisotropy is equal to that of a rock with the same porosity but made up of spherical pores. This is observed for fractures of both aspect ratios investigated here.

4. Discussion

4.1. P-Wave Velocity: EWAVE vs. MTEX-DEM

Pore-free rock P-wave velocities and anisotropy were calculated using the wave propagation model (EWAVE) and the Hill average effective medium model (MTEX), for samples cut perpendicular to foliation. The pore-free EWAVE and Hill average MTEX modeled velocities (velocities at 0% porosity in Figure 8A) are consistent between both methods, which is similar to observations by Zhong et al. (2014). As mentioned above, the non-indexed phases by EBSD are assumed to be mica. Based on Toy et al. (2008) and Dempsey et al. (2011), we also assume that most mica basal planes are parallel to foliation (80:20 parallel:random to foliation). Nonetheless, future studies should perform analyses capable of quantifying mica CPO on these mylonites to conclude on the true extent of the influence of mineralogy on the Alpine Fault rock anisotropy.

EWAVE and MTEX-GassDEM models are used to calculate the fractured rock P-wave velocity. A key difference between the model is that the GassDEM models a 3D pore, while the EWAVE model uses either natural rock fractures (e.g., extracted from micro-imaging techniques) or simplified 2D fractures, as done in this study. The simplified 2D fracture shape misses one dimension of the aspect ratio. For example, fractures with a 5:2:1 and 5:5:1 ratio would be simplified to a 5:1 2D ratio.

Individually these two types of 3D fractures (5:2:1 and 5:5:1) result in different GassDEM wave velocities, which also differ even more from the EWAVE modeled velocities using a 5:1 2D pore shape (Figure 10). The effect of pore aspect ratio on the fast velocities is insignificant. However, the 5:1 aspect ratio slow velocity deviates from the 5:2:1 and 5:5:1 aspect ratio trends by 12% and 5.3% at 5% porosity, respectively. As a result, the 5:1 ratio yields almost double the anisotropy of the 5:2:1 ratio.

Figure 10. Comparison between (A) modeled P-wave velocities and (B) P-wave anisotropy of the three rock models with three different pore aspect ratios: 5:1 (2D), 5:2:1 (3D), and 5:5:1 (3D).

A discrepancy between the elastic modeling using 2D and 3D data and their relationship has been previously reported (e.g., Saxena and Mavko, 2016; Hooghvorst et al., 2020). Because the 2D EWAVE model assumes a plane wave perpendicular to the EBSD model (Zhong and Frehner, 2018), the pore or fracture length perpendicular to the EBSD 2D model is thus assumed to be infinitely long. This is in contrast with the finite shape of the fracture used in the 3D GassDEM model. In this study, a simplified 2D pore shape is determined by reducing the 3D aspect ratio X:Y:1 to X:1. Having an infinitely long Y-axis in such a way contradicts this pore dimension reduction assumption and such infinite fractures exaggerate the effect of pores on the P-wave velocity. This explains why the modeled slow velocities from the 2D EWAVE model are always lower than those of the 3D GassDEM model. A key conclusion is that both the 2D and 3D pore aspect ratios strongly affect the modeled P-wave velocity and anisotropy. While both EWAVE and GassDEM models provide similar velocity predictions for pore-free rocks, a well-constrained pore shape in 3D is needed to accurately model P-wave velocity and anisotropy in a porous medium.

4.2. Pore Aspect Ratio Assumption

This study makes three assumptions regarding the pore aspect ratio used in EWAVE and GassDEM velocity modeling. First, all pores in a rock are assumed to have the same shape based on the average aspect ratio extracted from the micro-CT data. The average aspect ratio used in our models (5:2:1) is different from that of other works. For example, Vasin et al. (2013) uses a pore aspect ratio of 100:100:1 and 50:50:1 to model shale elasticity. Microfractures in shales also form parallel to layering, in particular, between clay platelets. The pore aspect ratio used in GassDEM modeling is referred to as an effective pore aspect ratio (e.g., Hong-Bing et al., 2013). It is the aspect ratio that best predicts the elasticity from experimental/field measurements at a fixed given porosity. Given that rocks contain various pore shapes and orientations, assuming one aspect ratio is idealistic and would not match all the real pore shapes in a rock.

Figure 4 shows that a pore network is composed of several connected pores of different shapes, sizes and orientations. The pore aspect ratio used in this work (5:2:1) represents the average pore shape of the whole rock, not individual pores which have a wide range of shapes (Figure 6). Finding a representative pore aspect ratio is challenging, but predicting pore shapes and volumes from laboratory data (Adam et al., 2020; Simpson et al., 2020) can be aided by modeling constrained micro-CT analysis for variable shapes and volumes of pores.

Second, we assumed that all lithologies of the shear zone and metamorphic Alpine Fault rocks in this study have the same pore shape. This assumption results in the same velocity relationship with porosity for all lithologies. Figure 6 shows that a similar dominant pore shape is estimated from the micro-CT data for all rocks. However, the pore shape distributions are different, pointing at the fact that the shape of individual pores could be different. In particular, the mylonite has flatness and elongation distributions that slightly differ from the protomylonite and ultramylonite.

Third, while modeling aligned fractures, we assumed that all pores have their longest axis aligned parallel to the rock foliation. Although this is the dominant direction for many microfractures, it is not realistic because not all pores align in the same direction in real rocks. A previous study on serpentinite (Kim et al., 2019) shows that a fracture effect on P-wave velocity is directionally dependent on fracture alignment. In our study, the fracture effect on velocity reduction is minimized when fractures align with wave propagation direction, and maximized when fractures align normal to the propagation direction. This results in the slow velocity being more susceptible to a change on fracture porosity than the fast velocity.

4.3. Velocity Comparison With Previous Studies

The pore-free modeled velocities of this study are compared to previous Alpine Fault EBSD-based velocity modeling works. The fast and slow pore-free velocities from this study (6.4 km/s, 5.5–5.9 km/s) are slightly higher than the velocities (6.2 km/s, 5.5 km/s) reported by Dempsey et al. (2011) but lower than the velocities (6.8–7.0 km/s, 5.7 km/s) reported by Adam et al. (2020). We suspect two possible causes for this discrepancy: mica volume and orientation. The rock samples Dempsey et al. (2011) analyzed contain 15-20 vol.% mica and their mica orientation is extracted from the EBSD data. Adam et al. (2020) worked with samples with a mica content of 35–40 vol.% and the orientation is assumed to align parallel to the rock foliation. In our study however, the samples contain 20–30 vol.% mica and only 80% of micas are assumed to align with foliation. This shows that understanding and quantifying volumes and orientations of mica-dominated rocks are critical to the estimation of seismic wave anisotropy.

We now compare our modeled velocity with the velocity measurements of Alpine Fault rock samples reported in previous studies (Okaya et al., 1995; Allen et al., 2017; Simpson et al., 2019; Adam et al., 2020). The velocities were measured under high confining pressure, thus assumed to be pore-free. The reported experimental fast and slow velocities range between 6.2–6.7 and 5.0–5.7 km/s, respectively. The velocities from this study fall within these ranges. Okaya et al. (1995) reports an increase in P-wave anisotropy toward the Alpine Fault plane. They interpret this increase in anisotropy as a result of increasing mica alignment due to metamorphism. Although we were unable to image mica orientations, we note that a 12% difference in mica volume fraction with 80% of micas being aligned with foliation can result in almost 5% increase in P-wave anisotropy (Figure 8D, Table 2). We further show that even small volumes of microfractures can significantly affect anisotropy. Adam et al. (2020) show that such fractures can remain open at 4 km depths at the Alpine Fault.

Field seismic data show that the rock in the vicinity of the Alpine Fault has a fast velocity of 6.2-6.5 km/s and a slow velocity of 5.6–5.7 km/s (Godfrey et al., 2002; Stern et al., 2007), which agrees with the velocities in this study. Additionally, the fluid-interpreted low velocity zone of 5.5 km/s velocity is also reported (Stern and McBride, 1998; Stern et al., 2007). Our porosity and fracture orientation models (Figure 9) predict the field velocities in the following sense: the fracture aspect ratio is 5:2:1 and 80% of fractures align with foliation while the rest are randomly oriented. It is important to recall that GassDEM models are non-unique, and different combinations of porosity, pore aspect ratio and orientation can result in identical velocity predictions.

5. Conclusion

In this study, P-wave elastic wave speeds in Alpine Fault rocks are modeled with and without microfractures. The four lithologies represent the protolith and the shear zone of the Alpine Fault: schist, protomylonite, mylonite and ultramylonite. Although samples were cut at 30^o and 90^o angles to foliation, the tiny mica minerals were not indexed in the EBSD analysis for both orientations. Through numerical modeling and laboratory data, the orientation of the micas is assumed to be 80% aligned to foliation and 20% random. Because all lithologies have similar mica volumes, the pore-free wave velocities and anisotropies are similar among the samples. However, we were not able to quantify mica CPO because of limitations of EBSD imaging of such small crystals. Therefore, it remains inconclusive whether variable mica CPO between the Alpine Fault lithology can have any significant influence on P-wave anisotropy.

The effect of fractures on the P-wave anisotropy of fault zone mylonite is the most important result of our study. Three sets of micro-CT data of protomylonite, mylonite and ultramylonite are analyzed to quantitatively extract porosity and pore shape information with a resulting average pore aspect ratio of 5:2:1 across all lithologies. Pore-free wave velocities are similarly predicted by EWAVE and MTEX. However, significant differences arise in the slow P-wave velocity prediction (i.e., anisotropy) when fractures are included as 2D rectangles (EWAVE) or 3D ellipsoids (GassDEM). The type of method used to model fractures can double the P-wave anisotropy in these rocks for a constant porosity. It is thus critical to understand how each numerical model approaches microfractures and to constrain shapes and volumes of such fractures with quantitative microstructural analysis.

Data Availability Statement

The original contributions presented in the study are included in the article/Supplementary Material, further inquiries can be directed to the corresponding author/s.

Author Contributions

JC prepared the thin sections, acquired EBSD data and performed the numerical models for this manuscript. LA developed the idea, co-wrote the manuscript and advised JC. VT and MO helped with the design and acquisition of the EBDS data and thin section preparation and revised the manuscript. BS and VT provided the micro-CT data. JS provided a modified GassDEM code and helped via discussions on numerical modeling. XZ provided advice on the EWAVE code and ways to modify it to include fractures.

Funding

We thank the Royal Society of New Zealand, Marsden Contract #14-UOA-028. The Avizo software and workstation employed were supported by Nvidia Corporation who donated the Titan X Pascal GPU, the Royal Society of New Zealand Rutherford Discovery Fellowship (RDF) 16-UOO-1602 and a subcontract to GNS Science (GNS-MBIE00056). Synchrotron and electron microscopic data acquisition was supported by RDF 16-UOO-1602, and SPring8 Proposals No. 2017B1378 and 2018A1506.

Conflict of Interest

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

Acknowledgments

We would like to thank Marianne Negrini for helping with the EBSD experiments. We also thank Andres Arcila-Rivera and Brent Pooley for helping with the sample preparation. Analyzing of the micro-CT data was assisted by Francesco Cappuccio.

Supplementary Material

The Supplementary Material for this article can be found online at: https://www.frontiersin.org/articles/10.3389/feart.2021.645532/full#supplementary-material

References