Measurement of 3D-Shape Preferred Orientation (SPO) Using Synchrotron μ-CT: Applications for Estimation of Fault Motion Sense in a Fault Gouge

We propose a 3D-shape preferred orientation (SPO) measurement method of rigid grains using synchrotron micro-computational tomography (μ-CT). The method includes oriented sampling, 3D μ-CT imaging, image filtering, ellipsoid fitting, and SPO measurement. After CT imaging, all processes are computerized, and the directions of thousands of rigid grains in 3D-space can be automatically measured. This method is optimized for estimating the orientation of the silt-sized rigid grains in fault gouge, which indicates P-shear direction in a fault system. This allows us to successfully deduce fault motion sense and quantify fault movement. Because this method requires a small amount of sample, it can be applied as an alternative to study fault systems, where the shear sense indicators are not distinct in the outcrop and the fault gouge is poorly developed. We applied the newly developed 3D-SPO method for a fault system in the Yangsan fault, one of the major faults in the southeastern Korean Peninsula, and observed the P-shear direction successfully.


Introduction
The presence of fault planes under the ground can be used as the basis to study the stress field at the time of paleo-fault activation and is one of the biggest factors affecting the stability of the ground evaluated in civil engineering and geotechnical engineering [1]. In addition, a development step of the fault is essential to understanding the tectonic revolution history. Accordingly, various methods have been proposed to utilize the micro texture of the rock to understand the major deformation mechanism of the fault. For instance, Lattice Preferred Orientation (LPO) analysis of clay minerals using X-rays and Anisotropy of Magnetic Susceptibility (AMS) analysis using magnetic properties of minerals are effective methods to understanding fault activation history [2][3][4][5]. In particular, recent studies that measured and analyzed AMS of a matrix composed of clay minerals included in the fault zone report that AMS is mostly affected by magnetic minerals such as hematite, but clay minerals are aligned similarly to the fault slip plane [6]. However, their object minerals were mostly limited to clay minerals in matrix, and there are only a few results of the texture analysis of minerals for fault gouges that exist in the brittle deformation zone.
The fault gouge is a kind of fault rock that has undergone an incohessive brittle deformation process, and it has visible fragments of less than 30% [7]. Therefore, rigid fragments and ductile clay matrix surrounding it are distributed in the fault gouge. Studies on behavior of rigid particles, To solve the aforementioned problems, micro-computational tomography (CT) imaging has been applied to various geology field [15][16][17][18]. The CT method can preserve the original sample, and a wide range of sample sizes can be measured. Recently, a research has been attempted to measure the shape of quartz, calcite, and pores in a strained shale through laboratory X-ray CT, and has been conducted on the strain difference given to rocks around fault zone [19]. However, unlike laboratory X-ray CT, when using a 3D imaging beam line installed in synchrotron, the energy level can be freely adjusted within a desired range to obtain a 3D image with better resolution. In this study, we attempted to obtain a high-resolution 3D µ-CT image of the fault gouge using a synchrotron accelerator by adjusting various factors such as wavelength and beam size. After the image filtering process, we attempted to automatize SPO measurement by coding the ellipsoid-fitting process of rigid grains and estimated the 3D fault motion sense quantitatively. Because the newly developed 3D-SPO method can automatically measure multiple rigid grains, this method can dramatically save time and improve the reliability of SPO results.

SPO Method Development
The newly developed 3D-SPO method comprises six steps-oriented sampling, 3D µ-CT imaging using the synchrotron beam line, image filtering and binarization, ellipsoid fitting, and SPO measurement. A flow chart of the steps is presented in Figure 2, and the details of each step are as follows.

µ-CT Imaging
We obtained µ-CT images at the beamline 6C Biomedical Imaging of the Pohang Light Source-II [20]. The beamline uses a multi-pole wiggler as its photon source and provides a monochromatic X-ray beam between 10 keV and 50 keV. For this experiment, we produced a 22-keV X-ray beam illumination using a double crystal monochromator based on silicon (111) reflection (DCM-V2; Vactron, Daegu, Korea). The specimen was mounted on an air-bearing rotation stage (ABRS-150MP-M-AS, Aerotech, Pittsburgh) 36 m from the source (Figure 3b). Behind the specimen, at a distance of 100 mm, an X-ray microscope was placed (Optique Peter, Lentilly, France) that comprised a 100-µm-thick CdWO4 scintillator facing the (010) direction (Miracrys LLC, Nizhny Novgorod, Russia), an objective lens of 4× magnification (UPLSAPO4X, Olympus, Tokyo, Japan), and a scientific CMOS camera (Zyla 5.5, Andor, Belfast, UK). The field of view (FOV) was 2.1 mm × 1.75 mm. For CT imaging, projections were acquired every 1.0 • till a 180 • rotation was complete. The exposure time was 0.5 s per projection.

Image Processing and Ellipsoid Fitting
Images acquired for each rotated angle during CT imaging were converted into a stacked 2D-image form through reconstruction (Figure 3b). For image preprocessing, CT reconstruction was performed using Octopus 8.9 software (XRE, Gent, Belgium) that supported filtered-back projection (FBP) reconstruction. The 3D image, which was converted into a stacking type of cube-shaped pixels, had a size of 2238 × 2238 × 1865, and the size of each pixel was approximately 0.94 µm. The resulting 3D images were processed and analyzed using Amira 6.5 (Thermo Fisher Scientific, Hillsboro, OR, USA). Before the image filtering process, each 3D-CT image was cropped in a cubic shape (800 µm × 800 µm × 800 µm) at the center of each sample, which had low noise ( Figure 4a). Several filtering processes were applied to allow easy differentiation between the grains and the matrix (or pore) in stacked 3D images. Among the various filters, smoothing and denoising filters were applied to remove noise, and a boundary sharpening filter was applied to clarify the boundary of the grains and the matrix (Figure 4b,c). In many cases, the filtering process was applied more than once, and the number of filtering times and type were performed depended on the type of sample (fault gouge, igneous rock, unconsolidated sedimentary rock, etc.). A non-local mean filter was used for denoising and smoothing effects and an unsharp masking filter was used for sharpening the boundary between the grains and the matrix (Figure 4c). In non-local mean filtering, small noises were removed, and the image was smoothened [21]. An unsharp masking filter is used to reclassify grain and matrix boundaries that are not easily distinguished because of excessive smoothing [22]. During the filtering process, the intensity range of the image is compressed compared to before processing. Finally, the intensity of the sample usually shows a distribution between 20 and 50 k. To distinguish by grain size, the grains in the image that are less than 100, 300, 500, and 1000 px in volume were removed from image.
The filtered images can be divided into grain, matrix, and pore, which are binarized in the next step. Binarizing the images involves visually comparing the original image and the binarized image and checking the grains by entering the appropriate boundary values (Figure 4d,e). Generally, among all components, grains have the highest intensity value and pores have the lowest intensity value. In this study, because SPO had to be measured by classifying only grains, image intensity was set by binarization by setting the intensity boundary value to classify grains. The standard intensity value that distinguishes the grains varied slightly depending on image. Therefore, it was better to visually compare the original image for each image (Figure 4d). This is because CT images measured by synchrotron show similar transmittance even though the mineral composition of each sample or the sample size is different. Therefore, it is better to set a different intensity for each image than to set one standard intensity for all samples. Additionally, small particles (<100 px) are not well processed with ellipsoid fitting or are difficult to separate even if noise is not deleted or grains are very small, which are problems faced during binarizing small particles ( Figure 4e). When measurement was performed on all samples including 100 px or less, the SPO directions of most small ellipsoids were in the directions of the principal X-, Y-, and Z-axes. This is because tiny objects with volumes of several tens of pixels are changed to simple ellipsoids with major axes in the direction of the X-, Y-, and Z-axes. Therefore, all grains with a volume less than the specific size (100 px in this study) in the binarized image were deleted. In addition, cut grains in the edge of the image were ignored in the analysis.
To estimate the orientation of each grain, the ellipsoid-fitting method developed which was adopted and be appropriate to fit an irregularly shaped object into ellipsoid [23]. By measuring the length of the x-, y-, and z-axes and the directions of the well-matched ellipsoid, the information of the aspect ratio (c/a), anisotropy (1-a/c), elongation (b/c), flatness (a/b), and the orientation of each axis can be obtained ( Figure 5). Fitting was performed for the binarized images using MATLAB code modified from the code of Kang et al. [23]. First, the code was developed assuming an ellipsoid with a volume of the object and to determine the length (a, b, c) of the principle axes (x, y, z) with the vector component and their orientations. Therefore, a covariance matrix (1) with a second moment µ value (2) was constructed. After the fitting is automatically applied, each grain is converted into an ellipsoid, from which the length and the direction can be obtained by measuring the length (a, b, c) of principle axes (x, y, z) and the orientation of the longest axis (z-axis).
x, y, and z are the axes of ellipsoid with vector components, and a, b, and c are the lengths of each axis.
V µ is the unit vector of the principal axes, λ 1 , λ 2 , and λ 3 are moments of inertia along the principal axes, and ν a , ν b , and ν c are matrix elements that represent the directivity of ellipsoids orthogonal to one another.
In addition, assuming the volume factor S as shown in Equation (3) below, the length of each major axis is calculated to be a = S·λ 1 , b = S·λ 2 , and c = S·λ 3 [23].
By coding the equations, we can make the best-fit ellipsoid for each rigid grain and can measure the length (a, b, c) and the orientation (x, y, z-axis direction, V a , V b , V c ) for the ellipsoid automatically.

SPO Measurement
The data obtained from each grain through ellipsoid fitting are the length (a, b, c) and the orientation (especially explained about V c in Figure 1b). We rearranged the whole data by the size (c), volume, and ellipsoid volume(4πabc/3), and the aspect ratio (c/a), which are crucial to interpret SPO results. The orientation (V c ) of the grain can be expressed by the direction vectors (Vx c , Vy c , Vz c ), which can be determined by measuring the trend (T), which is the rotated angle to the X axis in the XY plane, and the plunge (P), which is the angle between V c and Z-axis (Figure 1b). In practice, T and P values are apparent values (T a , P a ) measured with respect to reference plane. Therefore, by rotating these values using the orientation data of the sample, the real V c direction (T o , P o ) can be obtained. The T o and P o values of the target grains are plotted on an equal-area stereo-net with 5-degree intervals using MATLAB to present the distribution. The virtual measurement result of the SPO is shown in Figure 6a as an example, which shows how to display the fault plane (45 → 130) and the SPO (110 → 30) grains on the 2D stereo-net plane. In the same way, the ellipsoid's middle axis (b, V b ) and minor axis (a, V a ) were also investigated and plotted on a stereo-net. The internal angle (θ) between the tomographic plane and the SPO direction can also be obtained through the dot product of the vector. To consider the factors affecting the results, the overall SPO data were plotted on stereo-net after dividing into given size fraction (<5, 5-10, 10-15, >15 µm), given volume fraction (<100, 100-1000, 1000-10,000, >10,000 µm 3 ), and given aspect ratio (<1.5, 1.5-2, 2-3, >3). All data sorting processes were automated by MATLAB-based coding.

SPO Distribution Index
The distribution densities of directions of SPO values were presented as a contour map image. There are several indexes to express the distribution of 3D orientation such as Orientation Distribution Function (ODF), multiples of random distribution (m.r.d.), and mean tensor [24][25][26]. However, in the case of ODF containing a volume factor, there is a problem in that the weight is differently calculated for each grain. In addition, m.r.d. can express the relative distribution well through normalizing, but only presents quantitative values for densest points. Here we suggest an SPO distribution index (SDI) that can express the distribution density of SPO results quantitatively. The SDI indicates how the SPO data are scattered from the reference point, which is the most frequent direction. SDI is simply expressed as follows: SDI γ represents the maximum value of the angle between the γ% of grains surveyed and the most dense point. ϕ has a value between 0 and 180 • and is the angle between the direction which has the largest frequency (the densest point) and the direction of the major axis of the target grain. X ϕ is the number of objects in the direction of the major axis placed under the spatial angle of ϕ, and X T is the total number of objects. Therefore, the value of SDI ranges from 0 to 180 • , and it determines the distribution density, i.e., how many datapoints among the total are arranged in the corresponding direction. In addition, the distribution density of SPO can also be presented simply using SDI values with different γ values. The higher the SDI value for given same γ-values is, the more scattered the directions (V c ) of the z-axes of the grains are, whereas the lower the SDI is, the more concentrated the arrangement of grains in a specific direction is.

Sample Preparation
The gouge sample was collected from a outcrop of Yangsan fault zone, which is one of the largest fault zones in the southeastern parts of the Korean Peninsula and has the longest and most prominent linear structure in the NNE direction of strike [27][28][29]. The sampling site (BK) is located in the northern part of the Yangsan fault, which penetrates through the fault core with a typical fault gouge. The fault zone is composed of a core zone (9−10 m in thickness) and surrounding damage zones (>200 m in thickness) [30]. In the fault core zone, lens shaped rock fragments of several centimeters to one meter in size appear, indicating a dextral sense of shear. Most clay minerals contained in the matrix are illite and chlorite, and some smectite is also present [30]. Sim et al. determined the fault activation ages were 20 and 26-28 Ma using the Illite age analysis (IAA) [31]. The fault plane observed in the outcrop stands vertically, and the direction of the fault plane is measured as N45• E 85• SE. In the fault core zone, the sample was taken from the fault gouge, and it is blue-grey colored and less than 50 cm wide.
To obtain samples for the experiment, the iron core was inserted into a fault gouge using a hammer, and then the strike and dip direction of the inserted surface was measured (Figure 3a). Before the sample was taken out, the strike direction of the injected iron core was marked. The collected samples were sealed during transport from the field to the laboratory, where the samples were dried for more than 24 h. Cold-mounting was performed using glycol methacrylate (GMA) resin which has similar viscosity with water at room temperature, because it was difficult to maintain the shape of the sample near the unconsolidated state by making slabs or thin sections. A GMA resin was used to solidify the sample at 60 • C under vacuum. For some samples, the cold-mounting process was repeated several times to fabricate a complete solid form of the specimen with the resin penetrating deep into the pore cavity. After cold-mounting was completed, the sample was cut into a size of 1 mm × 1 mm × 2 mm.

SPO Data of Fault Gouge Sample
µ-CT images of the fault gouge were obtained using the synchrotron beam line of the Pohang Light Source-II [20]. After image filtering for the µ-CT image, 3D-SPO estimation was performed by measuring the length (a, b, c) and direction vector (V a , V b , V c ) for the ellipsoid-fitted grains of 7636 particles. The 3D-SPO results of major axis were plotted with red color on the stereo-net with 5-degree intervals. In addition, the results of measuring the direction of the middle axis and minor axis of the ellipsoid are shown together in blue and green on the stereo-net, respectively (Figure 6b). Examples of data are shown in Table 1. The grain size (c), which is the length of the z-axis of the fitted ellipsoid, ranged from 3.2 to 82.4 µm, and the grains of 85% have less than 10 µm sized. As the grain size increased, the number of grains tended to dramatically decrease. The range of the aspect ratio (c/a) is various, but the numbers 2-3 and >3 are the greatest. More detailed information of the investigated grains, such as anisotropy, flatness, elongation, graph of aspect ratio to grain size, and Scanning Electron Microscopy (SEM) images is given in Supplementary Materials (Figure S1-S3). SPO results were presented as a contour map on the stereo-net (Figures 6b and 7). The point of the highest SPO density distribution was found to be 028 → 24 and has 17.0 of SDI 20 values. Various SDI examples with γ-values varying from 10 to 50 with 10 intervals for the 3D-SPO results are presented in Supplementary Materials Table S1. The SDI 20 values are presented on the contour plot to show the degree of distribution of the SPO results ( Figure 6b). Finally, the angle (θ) between the fault plane and SPO direction was calculated from the dot product of the vector to verify whether it corresponds to the P-shear direction of the fault plane (Table 1).  The 3D-SPO results were classified into several groups by grain size, volume, and aspect ratio. the contour plots of each results are separately presented in Figure 7. The directions of the densest points for the each SPO data by grain size appeared in the order 028 → 22, 029 → 22, 023 → 20, and 025 → 22, with an SDI 20 value of almost 20 (Figure 7a-d). The directions were similar to the overall direction (028 → 22) as shown in Figure 6b, whereas the SDI 20 values were larger than the overall value (17.0). The directions of the densest points for the each SPO data by grain volume appeared in the order 030 → 24, 025 → 23, 026 → 23, and 022 → 20, with an SDI 20 value of 16 to 21 (Figure 7e-h). For the SPO data by each aspect ratio, the directions of the densest point are different and range from 018 → 14 to 030 → 29 (Figure 7i-l). These results were also slightly different; the overall direction is 028 → 24, and the SDI 20 values are approximately 20.

SPO Data Interpretation and Implication
In this study, we successfully obtained 3D-SPO results from synchrotron µ-CT images. This developed 3D-SPO method allowed us to investigate the hidden vertical displacement of the strike-slip fault and the unmeasurable horizontal displacement of the vertical fault in the millimeter and micrometer scales. The SPO results of the middle and minor axes appeared scattered on a plane perpendicular to the SPO results of the major axes. Relatively, the data in the middle axis showed a dense form at two points, and these results, along with the characteristics of grains included in Supplementary Materials, could be considered as further research topics. On the other hand, we tried to identify the factors controlling SPO results. The case study showed that the spatial distribution of SPO looks more dispersed as grain sizes increased (Figure 7a-d). In particular, for the grains of sizes of <5 µm (2834) and 5 to 10 µm (3652), the SPO results were radially well-distributed on the contour plot (Figure 7a,b), whereas the contour plots of grains >15 µm (327) were relatively irregular (Figure 7d) with no distinct changes in the SPO direction. This could have been caused by the abundance of the size fractions of the target fault gouge. Although it was analyzed as above in the data classified by grain size, when the sufficient SPO data were collected, the difference in the tendency of the direction in each grain size was not large. This can also be confirmed in the SPO data classified by volume. The total SPO data is well matched the case of 100-1000 µm 3 (5860) where most of the data is concentrated (Figure 7f). However, the other three cases, which have little data, have a larger SDI value or even some scattering points (Figure 7e). Therefore, the results of analysis of the contour plot and SDI distribution by particle size effect show that when studies are conducted on a small number of particles, the SPO data may be inaccurate. As a result, if multiple particles are included, regardless of particle size, the accuracy and precision of SPO results can be improved. As the SPO results were classified by aspect ratios, all results showed a similar direction for grains with similar SDI values, except for aspect ratios <1.5 (Figure 7i-l). The SDI 20 value of the grain fraction with <1.5 aspect ratio increased slightly. This suggests that if grains are shaped like a sphere, they show a low directional density of the major axis. This is because particles having a relatively low aspect ratio are less likely to be arranged in the P shear direction during rotation of the rigid particles. This is accordant with the previous research result about rigid body rotation of particles with different aspect ratios using computational modeling [32]. However, SPO data of grains appearing in the natural fault zone does not show a significant difference from the total data when it has an aspect ratio of 1.5 or more. Therefore, checking the measurement of grains with a relatively high aspect ratio results in better accuracy and precision of 3D-SPO results.

Estimation of Fault Motion Direction
According to previous research, the rigid particles, such as quartz, are free to rotate within a matrix of clay and other fine-grained material without interference from similar sized grains [13,14,32]. In addition, the SPO direction of the rigid grains distributed in the fault gouge reflects the movement of the fault plane and is in the P-shear direction [14]. Based on this concept, we can inversely estimate the movement direction by the shear of a given fault plane using the SPO result. In the 028 → 24 direction, the SPO density point of grains investigated in this case study showed an angle about 20 • with respect to the fault plane (N45 • E 85 • SE). Assuming that the developed fault plane moved in the form of a dextral slip, the result corresponded to the P-shear direction. Therefore, it can be confirmed that the fault plane has a fault movement in the form of dextral slip. Similarly, the SPO is distributed at low angles with respect to the vertical displacement, indicating that the hanging wall of the southeast was uplifted, and finally, it underwent reverse fault movement. Thus, the fault can be of the dextral and reverse types, with relatively prominent movement of horizontal displacement. This result is consistent with the movement determined by shear sense indicators observed in outcrops [30]. However, in the previous study, only a simple information of fault movement was defined. In contrast, the 3D-SPO results of this study determined the direction of the spatial fault motion in 3D.
In previous studies, the SPOs were measured only on 2D images of specific cross sections, and the number of investigated grains were not enough, thus limiting the estimation of the fault slip direction and decreasing the reliability of SPO data [14,33]. The SPO results using 2D images were not sufficient, and only half of the information for the faults with both vertical and horizontal displacements was obtained. In field surveys, the fault outcrop is a particular cross section, which provides a limited data for fault activation and can cause errors due to regional differences. In contrast, the newly developed 3D-SPO method using 3D µ-CT can overcome the various limitations in estimating fault motion directions and can be applied to investigate fault gouges with a thickness of less than 1 cm. Therefore, SPO measurements using µ-CT are the most versatile evidence for determining the fault motion direction in outcrops with fault gouges.

Conclusions
The developed 3D SPO method using synchrotron µ-CT imaging was applied to define the fault motion direction of rigid grains generated by fault activation in a brittle shear zone case. In this study, we confirmed that the dextral strike-slip-dominated fault movement with some reverse sense was defined. Moreover, the accuracy and precision of the data considerably increased by measuring a greater number of fine-silt fraction's SPO (<5 µm and 5-10 µm) with an aspect ratio of >1.5. In conclusion, the 3D-SPO method can quantitively restore the motion of fault and can overcome the limitations of 2D-SPO studies both at the field scale and at the laboratory scale with high reliability. In addition, this method is versatile enough to provide basic data on future studies on earthquakes for quaternary active faults with thin gouges.