Observation of Temporal Variations in Seismic Anisotropy Within an Active Fault‐Zone Revealed From the Taiwan Chelungpu‐Fault Drilling Project Borehole Seismic Array

Temporal fault‐zone observations are important to better understand the evolution of fault structure and stress configuration. However, long‐term monitoring in the fault zone is rare after a large earthquake. Here, we use seismic data in the fault zone at 1‐km depth from the Taiwan Chelungpu‐fault Drilling Project to study long‐term anisotropy after the 1999 Mw7.6 Chi‐Chi earthquake. The direct S‐wave splitting measurements resolve the overall weak anisotropy in the shallow crust. In order to resolve fault damage zone anisotropy, we perform coda cross‐correlation technique for 794 microearthquakes between 2007 and 2013. We estimate the temporal change in background shear‐wave velocity, fast shear‐wave polarization direction (FSP), and strength of anisotropy (Aani) in the fault damage zone. We show the average FSP direction is N93°E with a significant Aani of about 12%, likely due to the pervasive vertical microcracks created after the earthquake. Temporal variations of anisotropy exhibit seasonal variation with periodicity every 9–12 months that correlates with rainfall events. Furthermore, long‐term anisotropy shows a gradual rotation of FSP direction of about 15° during the first 4 years of observation. At the same time, the strength of anisotropy reduced from 17% to 10% and shear‐wave velocity increased, suggesting the fault healed after the earthquake. This study reports in‐situ evidence for two key observations: (a) long‐term, fault‐zone healing after a major earthquake and (b) modulation of 1‐km deep fault‐zone properties by surficial hydrologic processes. These observations may provide constraints on the response of the fault damage zone in the interseismic period.

Seismic anisotropy is one useful tool to observe the local stress field in the fault zone and its temporal evolution. Seismic anisotropy is a phenomenon in which seismic waves travel at different velocities as a function of wave orientation. This variation in seismic speed with direction can be due to intrinsic rock properties (e.g., Johnston & Christensen, 1995), and/or can be stress-induced (Boness & Zoback, 2004;Hung et al., 2009;Zatsepin & Crampin, 1997). Stress-induced anisotropy results from the differential stress which generates vertically-aligned fractures orienting to the horizontal maximum principal stress (S Hmax ) direction (Crampin, 1987(Crampin, , 1994. Therefore, when a shear-wave propagates in an anisotropic crust, the incoming wave splits into two quasi-shear waves traveling at different velocities (i.e., shear-wave splitting, SWS). Two observational factors are indicative of seismic anisotropy: (a) the fast shear-wave polarization (FSP) direction which is aligned with the predominant fracture orientation and (b) the strength of anisotropy which is recorded by the time delay between the fast and slow shear-waves. In general, the time delays are related to the magnitude of differential stress in the fault zone (Zatsepin & Crampin, 1997).
Several factors can cause anisotropy to vary in space and time, including fault fabric and fracture alignments Liu et al., 2008;Zhang & Schwartz, 1994), overpressured fluids (Angerer et al., 2002;Crampin et al., 2002;T. Li et al., 2019), coseismic stress variations (Saiga et al., 2003), and fault-zone damage (Tadokoro & Ando, 2002). Temporal variations in anisotropy are associated with coseismic damage in the fault zone from newly generated fractures (Durand et al., 2011;Tadokoro & Ando, 2002). Some of the newly generated fractures may heal in the postseismic period, leading to a recovery of fault-zone seismic velocity or anisotropy (Kaproth & Marone, 2014;Tadokoro & Ando, 2002). This healing process has been observed to be a logarithm function of time that can last for several years (Baisch & Bokelmann, 2001;Kaproth & Marone, 2014;Schaff & Beroza, 2004;. Tadokoro and Ando (2002) studied the postseismic fault-zone behavior after the 1995 Kobe, Japan, earthquake, and suggested that the temporal fault-zone evolution may be inferred from anisotropy changes. However, such observations of long-term seismic anisotropy in a fault-zone after a major earthquake are still rare. In this study, we take advantage of the seven-level vertical seismic array of Taiwan Chelungpu-fault Drilling Project Borehole Seismometers (TCDPBHS, hereafter) to examine the temporal seismic anisotropy after the 1999 Mw7.6 Chi-Chi earthquake. The TCDPBHS was installed across the fault-zone, where a primary slip zone was identified at the depth of 1,111 m from the retrieved fault core (Ma, 2021;Ma et al., 2006). This unique observational dataset allows us to answer the following questions: how does the Chelungpu fault evolve in response to postseismic fracture healing and background stresses (tectonic and pore pressure)? And how does subsequent healing continue in the interseismic phase? We perform temporal anisotropy analysis in the Chelungpu fault-zone to monitor its healing process. Specifically, we will primarily focus on the anisotropy analysis from coda cross-correlation to monitor the fault-zone anisotropy and its temporal variations. In addition, we analyze earthquake shear-wave splitting to gain insight into crustal anisotropy in the region.

Waveform Data and Microearthquakes
The TCDPBHS, consisting of seven vertical-leveled seismic stations (short-period velocity and natural frequency 4.5 Hz), was installed at depths of 946-1,274 m (from BHS1 to BHS7) with an interstation spacing of 50-60 m ( Figure 1a). The BHS4 is at a depth of 1,110 m, close to the primary slip zone of the Chi-Chi earthquake at 1,111 m. The stations BHS1 to BHS4 (950-1,110 m) were placed within the damage zone of the hanging wall of 10.1029/2021JB023050 3 of 19 the Chelungpu fault, mostly in the Cholan formation (fm.) and Chenshui Shale fm. The BHS5 to BHS7 stations were placed in the footwall, predominantly in the Chinshui Shale fm (Figures 1a and 1b, A. T. S. Song et al., 2007). The TCDPBHS has been functioning since late-2006, providing continuous waveform data. The sampling rate of the data was initially set to 1,000 Hz from November 2006 to December 2007, and it was decimated to 200 Hz in January 2008. To unify the data analysis, we down-sampled the data in 2007 to 200 Hz for consistency. All the data are corrected for the instrumental response, Galperin angle, and the sensor orientation (see details in Lin et al., 2012).
In this study, we focus our analysis on data from 2007 to 2013. Due to unlock of the GPS clock from mid-2010 to 2011, accurate traveltime information was not attainable during that time period. Therefore, we do not apply further analysis against data in this period. Among the seven stations placed in the damage zone, BHS1 and BHS4 are the most reliable. Data from other stations is contaminated by instrumental noises (Y. J. Wang et al., 2012). As a result, we use waveforms recorded at BHS1 and BHS4 to analyze the in-situ fault zone anisotropy.
As shown in Lin et al. (2012), regional microearthquakes are mostly distributed at around 10 km depth in a confined area (Figure 1c). This consistent source region provides a stable azimuthal angle to monitor the longterm evolution of seismic structure and to infer the state of stress in the fault system. To investigate the temporal evolution of seismic anisotropy, we first extend the existing microearthquake database using the semi-auto picking algorithm  to detect possible events from continuous BHS1 and BHS4 waveforms. We manually checked the signals of identified microearthquakes to verify the detections. Subsequently, we obtain the location of these microearthquakes by analyzing the shear wave polarization described in Supporting information and Figure S1. The magnitude of detected microearthquakes, M we , is determined by the empirical relationship  . The brown dashed line indicates the primary slip zone at a depth of 1111m . (b) cross-section of the drilling site and the geometry of the Chelungpu thrust fault dipping (∼30°) to the east. This figure is modified from Hung et al. (2009). (c) A local map shows the TCDP drilling site (black triangle) and the weather station C0F97 (black square) for rainfall data collection. Local microearthquakes from November 2006 to December 2007 identified by Lin et al. (2012) are shown in solid dots, whereas crosses are the microearthquakes used for coda cross-correlation between 2007 and 2013. Leftbottom subfigure indicates our study area's location relative to a regional map of Taiwan.
For the anisotropy analysis, we setup three criteria for the event selection: (a) events with S-P travel time less than 2 s (hypocenter distance less than 15 km) to reduce the complex path effect; (b) events with the M we greater than 0.8 to ensure high-quality waveforms; (c) events with the shear-wave incident angle less than 40° to avoid S-to-P conversion (Crampin, 1989). Based on the considerations given above, 794 out of 21,014 microearthquakes from 2007 to 2013 are used for the fault-zone anisotropy analysis. Figure 2 shows the selected events and their temporal and spatial distributions with respect to the TCDP site. Seismicity is concentrated in the northeastern area, generally consistent with the observations shown in Lin et al. (2012) (Figures 1c and S1).

Estimation of Seismic Anisotropy
Seismic anisotropy may be inferred from earthquake-generating SWS which is useful to study crustal structures and stress configurations (e.g., Chang et al., 2009;Cochran et al., 2003;Crampin et al., 2002;Kuo et al., 1994;Liu et al., 2008;Rau et al., 2000;Savage et al., 1990;Zhang & Schwartz, 1994). More recently, anisotropy analyses are carried out using ambient noise or coda-wave interferometry (e.g., L. W. Chen et al., 2017;Lewis & Gerstoft, 2012;Miyazawa et al., 2008;Nakata & Snieder, 2012;Tonegawa et al., 2013), which retrieves the interstation Green's function from cross-correlation of noise/coda waveforms. Unlike direct waves which only sample the medium once, coda waves are generated by multiple reverberations from scattering. The multiple samples generated by coda waves may enhance the resolution of the physical structure (Grêt et al., 2006). Such a technique provides a greater capability in constraining seismic anisotropy at the local scale (e.g., ∼100 m, L. W. Chen et al. (2017); Lewis & Gerstoft (2012)). In the following sections, we address the methodology for anisotropy estimation via coda cross-correlation and direct shear-wave splitting measurements, as well as temporal monitoring.

Cross-Correlation From Coda Waves
The cross-correlation of seismic noise or coda waves between two stations has been shown to be a feasible way to retrieve the interstation Green's function, that is, 1,2 ≈ − 1,2 ( 1, 2, − ) + 2,1 ( 1, 2, ) , where 1,2 is the cross-correlation function of two displacement waveforms, and ( 1, 2, ) represents the Green's function (Compillo & Paul, 2003;Sabra et al., 2005;Snieder, 2004). This technique has been widely used to investigate crustal structures (e.g., L. W. Chen et al., 2017;Lewis & Gerstoft, 2012;Nakata & Snieder, 2012;Tonegawa et al., 2013). However, this method assumes that noise sources are evenly distributed in space. In reality, the distribution of noise sources may vary over seasons (e.g., Landès et al., 2010), resulting in some degrees of variation in the cross-correlation function, and leading to inaccurate travel time measurements (Froment et al., 2010;Tsai, 2009). On the other hand, coda waves are generated from multiple scattering (Aki, 1969), producing a more diffusive wavefield and a more robust retrieval of Green's function (Froment et al., 2010). In order to mitigate the potential bias in phase traveltime due to uneven noise source distribution, the coda wave cross-correlation is used in this study. The spatial distribution of selected earthquakes with an incidence less than 40° is generally consistent (i.e., in the northeastern TCDP, Figure 2) overtime, and permits a robust Green's function construction (e.g., Emoto et al., 2015).
We define a 2-s coda window starting 0.2 s after the S arrival from the microearthquakes, where most of the coda perturbations remain strong ( Figure S2). We adopted the program MSNoise (Lecocq et al., 2014) to compute the coda cross-correlation function (CCCF) of the horizontal component between BHS1 and BHS4, with BHS1 being the virtual source and BHS4 the receiver. The maximum time lag of the CCCF is set at ±1s. The preprocessing whitening is applied in the frequency domain.
To ensure the equal contribution of individual CCCF to the stacked CCCF, we use the phase-normalized stacking method, where the amplitude of each CCCF is normalized by its maximum amplitude in the frequency band of 5-30 Hz before stacking. The phase arrival time can be determined from the peaks of the stacked CCCF. We test the stability of phase arrival against different numbers of stacking and find that the stacking number of 40 permits a stable CCCF phase arrival (hereafter simple stacking, Figure S3a). The phase arrival time between BHS1 and BHS4, separated by 164 m, is about 0.08 s, corresponding to the shear-wave velocity of about 2 km/s ( Figure S3b and S3c). This velocity estimate is also similar to the result from ambient noise cross-correlation at TCDPBHS (Hillers et al., 2012), showing agreement with the well-logging measurement (H. Y. Wu et al., 2007, Figure 1a). These consistencies suggest that our stacked CCCFs represent the shear-wave Green's function between BHS1 and BHS4.
To determine seismic anisotropy between BHS1 and BHS4, we repeat the above calculation by rotating the horizontal components from 0° to 180° with a 5° increment to measure azimuthal variations in the S phase arrival time. To estimate the anisotropy parameters, we apply a cosine fitting defined as where T (θ) is the travel time measured along the azimuth of the horizontal component θ, T 0 is the isotropic shearwave travel time, T 1 is the anisotropic coefficient and is the FSP (Alfold, 1986;Nakata & Snieder, 2012). The strength of anisotropy A ani is defined as A ani = (T slow -T fast )/T 0 × 100%, where T slow and T fast are shear-wave travel time measured in the slow and fast polarization directions, respectively. An example result from stacking 40 CCCFs is shown in Figures 3a and 3b, where the FSP is N105°E and the A ani is ∼13.9%, corresponding to T slow -T fast ∼ 0.01 s. The time window for the stacked CCCF in this example covers data from about 5 months (4 January -28 April 2007). To evaluate the robustness of A ani and , we repeat the CCCF calculation against interpolated data (spline line interpolation with a sampling rate of 10,000 Hz). We found that the results estimated based upon CCCFs from non-interpolated data and those from interpolated data are very similar ( Figure S4).
Furthermore, we use the bootstrap method (Efron & Tibshirani, 1994) to verify the reliability of anisotropy parameters. Specifically, we randomly select 40 CCCFs from all CCCFs within a given time window to generate a new group for stacking. In each group, CCCFs are allowed to be selected repeatedly. Such a process is repeated 200 times, and we calculated the mean and the standard deviation of T 0 , T 1 , and to understand the uncertainty of the estimated anisotropy. An example after bootstrap method from the same time window (4 January -28 April 2007) is shown in Figure 3c, where the uncertainties of the FSP direction and A ani are about 5° and 5.7%, respectively. The estimated anisotropy parameters are very similar to those determined by the simple stacking method (Figure 3b), suggesting the reliability of these estimates with variances.

Measuring Shear-Wave Splitting (SWS)
In addition to the coda cross-correlation method discussed in Section 3.1, we also conduct SWS measurements against the direct shear-wave from microearthquakes. Here we adopt the cross-correlation method to search for the optimal anisotropy parameters (Bowman & Ando, 1987). First, the S waveforms are band-pass filtered between 1 and 20 Hz to minimize noise contamination. The horizontal-component seismograms are rotated between 0 and 180° in a 5° interval before cross-correlation. The maximum delay is set to 0.05 s, but it can be up to 0.1 s in some cases. We apply the grid-search approach to identify the waveform pair with the greatest normalized cross-correlation coefficient. The normalized cross-correlation value (CCV) is defined as where x i and y i represent S waveforms at two horizontal components, and n is the cross-correlation window length. We found that filtered S waveforms with one or two wave cycles are sufficient to produce a robust measurement.
We use the Fisher Transformation (e.g., Kreyszig, 1970) to estimate the uncertainty of CCV, assuming the two variables, rotation angle and time delay, follow the bivariate normal distribution. We define 95% confidence intervals and contour this interval in the cross-correlation domain to circle the range of uncertainty (see details in Rau et al., 2000). The error bar of the FSP direction and time delay is determined by taking the difference between the greatest uncertainty and the optimum parameters within this range. Figure 4 illustrates an example of SWS measurement. The FSP directions at BHS1 and BHS4 show similar orientations of about 115-120° and an identical time delay of 0.045 s.

Measuring the Temporal Anisotropy
We use moving window stacking of the CCCFs to estimate temporal anisotropy parameters. The moving window contains 40 CCCFs, and it is shifted by one CCCF in each step. Each moving window spans roughly 4 months (i.e., 80-200 days with a mean of 118.4 days). We stack over CCCFs in different azimuths and apply Equation 1 to obtain the anisotropy parameters in each moving window. We discard results from time windows that cover time gaps over 100 days (e.g., due to GPS unlock, Figure 2), to preserve a smooth anisotropy pattern. Since the data sampling rate was 200 Hz, A ani smaller than 6.25% (or time delay = 0.005s) is not included in the subsequent analysis. Consequently, the results from 715 moving windows between 2007 and 2013 ( Figure 5) are used to elaborate on temporal change in anisotropy. For the seismic anisotropy in the bulk upper crust from SWS measurements, the moving window technique is applied to the temporal anisotropy as well. We set up a shorter moving-window length to analyze the temporal behavior of anisotropy because the SWS method resolves anisotropy parameters without a stacking process. We use moving windows that include 10 events (mean of about 36 days) to calculate the average anisotropy parameter FSP and time delay dt in each window. The average FSP direction (FSP avg ) is calculated through the following equation where the is FSP of ith earthquake, and k is the total number of earthquakes. The average time delay is simply ∑ . The moving window is shifted by one event in each step. Equation 3 refers to the circular mean, and is useful in presenting the overall anisotropy in a given area (e.g., Gerst & Savage, 2004).

Results and Discussion
The methods mentioned above allow us to study the temporal anisotropy at the crustal-scale and the fault damage zone scale. In the following sections, we elaborate on our findings. Specifically, we will start with the anisotropy from SWS results, and provide a general picture of the anisotropy in the bulk upper crust. Then, we will focus on the interpretations of the temporal anisotropy in the fault damage zone based on the CCCF analysis.

Seismic Anisotropy in the Bulk Upper Crust From SWS
We find the FSP direction from SWS ranges between 90 and 120° for BHS1 and 90-160° for BHS4, respectively, with very similar time delays of about 0.02-0.04 s ( Figure S5 and Figure 6). The SWS analysis cannot resolve fault-zone anisotropy between BHS1 and BHS4 (i.e., ∼164 m), because the dominant frequency of the S wave in SWS measurements is ∼9 Hz and the wavelength is on the order of 200 m. Furthermore, seismic scattering associated with the fault zone could also make it challenging to obtain a stable splitting time measurement (e.g., Aster et al., 1990). Instead, we analyze SWS data from 2007 to 2008 to understand the overall background anisotropy in the shallow crust. SWS analysis suggests the strength of anisotropy in the shallow crust (about 10 km) beneath the Chelungpu fault is consistently weak over the time period between 2007 and 2008 (∼0.025s, blue lines in Figure S5). For comparison, Chang et al. (2009) found that the majority of SWS time delays in the Western Foothill and mountain areas of Taiwan are about 0.05 s over 20 km in the upper crust. The obtained time delay in our study generally agrees with that shown in Chang et al. (2009). The FSP avg direction over the entire dataset is about 101 ± 30° at BHS1 and 107 ± 42° at BHS4, respectively. These estimates are in general compatible with the tectonic convergence direction (Yu et al., 1997).  diverse incident angles and azimuths. As the FSP and time delay depending on the ray azimuth and incidence (Crampin, 1985), it is not surprising that FSP directions from these two methods are not exactly the same.
The time delay from SWS is ∼0.025 s averaged over the 9-14 km between the microearthquake source region and TCDP seismometers. By contrast, the time delay from CCCF analysis is ∼0.011 s averaged over the 164 m between BHS1 and BHS4. The time delay measured inside the damage zone indicates strong anisotropy over such a short distance. This is a significant difference in the strength of anisotropy as compared to the weak anisotropy in the background crust. Chang et al. (2009) interpreted the weak crustal anisotropy to be the result of overlays of diverse fossil fabrics and fractures in past orogenic processes. On the other hand, the strong anisotropy observed in the damage zone by CCCF analysis could be the result of pervasive fractures with similar alignments documented by in-situ well logs from the Chelungpu fault damage zone (Hung et al., 2009;H. Y. Wu et al., 2007; see Section 4.2). The CCCF results between BHS1 and BHS4 constrain in-situ anisotropy within the fault damage zone, whereas SWS results are more indicative of the bulk seismic anisotropy in the upper crust.

Seismic Anisotropy in the Fault-Zone
Based on Section 4.1, we focus on the anisotropy in the fault damage zone from CCCF results, and monitor its temporal behavior. As shown in Figure 5, the FSP direction varies between 50 and 120°, with an average of The observed strength of anisotropy and averaged FSP direction from CCCF analysis could be the result of layering in sedimentary rocks, alignment of anisotropic minerals, and/or stress-induced microcrack orientation. Determining the primary mechanism for the observed seismic anisotropy is key to our understanding of whether or not the anisotropy is controlled by fractures. Between BHS1 and BHS4, Choloan fm. is dominated by thick sandstones with alternating layers of sandstone, mudstone, and siltstone, whereas the Chinshui Shale fm. is dominated by siltstone with thin, alternating sandstone-siltstone layers (Hung et al., 2009;Louis et al., 2008). As noted by Hung et al. (2009) and Y.H. Wu et al. (2008), the anisotropy due to alternating formation layers with relatively small dipping angles (i.e., <30-40°) is minimal. This is also consistent with the observations that the FSP direction from the DSI both deviate from the strike of the bedding plane by about 90° (Hung et al., 2009;Y.H. Wu et al., 2008). Therefore, sedimentary layering is likely not the primary cause of the observed anisotropy.
Stress-induced fracture orientations may play an important role in seismic anisotropy as the strength of anisotropy increases with increasing fracture density (Crampin, 1994). Near TCDP, H.Y. Wu et al. (2007) observed the fault-zone fractures by using in-situ well-logging methods and found 63% of the fractures are concentrated at depths between 500 and 1,100 m. In addition, Y.J. Wang et al. (2012) examined the attenuation factor in the fault zone and discovered a low Qs of about 21, potentially associated with a high fracture density. Using the core samples from TCDP, Louis et al. (2008) showed that the P-wave anisotropy in the Cholan fm., mostly sandstones, is likely related to the alignment of microcracks in the maximum stress direction. On the other hand, the P-wave anisotropy in siltstone samples, mostly from the Chinshui Shale fm., is dominated by the preferred orientation of phyllosilicate minerals, suggesting that the alignment of grains/minerals may also contribute to the observed anisotropy (e.g., Johnston & Christensen, 1995).
Despite the fact that the Chinshui Shale fm. is thicker than the Cholan fm. by a factor ∼2 between BHS1 and BHS4 (Figure 1a), the strength of anisotropy from the sandstone samples is at least 3 times higher than that from the siltstone samples (Humbert et al., 2012;Louis et al., 2008). Humbert et al. (2012) also showed that the S-wave anisotropy in the sandstone and siltstone samples are about 15% and 3.5%, respectively. It's worth noting that the observed FSP direction is compatible with the alignment of vertical cracks and the FSP direction in the sandstone samples (i.e., ∼N105°E, Humbert et al., 2012;Louis et al., 2008). In addition, the observed FSP direction is almost 90° from the FSP direction associated with the siltstone samples (i.e., fault-strike parallel, Humbert et al., 2012). While both the Cholan and Chinshui Shale fm. could contribute to the observations, it is expected that the Cholan fm., mostly sandstones with strong anisotropy, contribute more substantially to the observed anisotropy parameters than that from the Chinshui shale fm. (e.g., Silver & Savage, 1994), especially when the preferred orientation of minerals is weak (Y.H. Wu et al., 2008).

Temporal Variations of Seismic Anisotropy in the Fault-Zone Between 2007 and 2013
As discussed earlier, there is strong evidence to support a highly fractured Chelungpu fault-zone. Zatsepin and Crampin (1997) suggested that the horizontal differential stress (i.e., ∆σ = S Hmax − S hmin , where S hmin is the horizontal minimum principal stress) plays a critical role in stress-induced anisotropy. As the preferred orientation of phyllosilicate minerals is likely to be stable over the time scale of seismic cycles, the change of anisotropy may be predominantly associated with the change of ∆σ and the orientation of microcracks. Any stress change, either in direction or in magnitude, can modify the ∆σ and consequently cause the anisotropy to vary. In this instance, we interpret temporal variations of anisotropy in the fault zone to reflect fracture evolution.
Here, we first show three CCCFs as a function of calendar time against azimuths (Figures 8a and 8b). The waveform arrival time shows substantial variations as a function of wave polarization direction. Waves polarized in the northeastern-southwestern quadrant show the greatest arrival time variations, while waves polarized in the WNW-ESE quadrant show the least (Figure 8c). This is expected given that the core samples show vertical cracks are aligned along the WNW-ESE direction (i.e., N105°E, Louis et al., 2008). As the crack density increases, we expect that the shear-wave speed varies strongly with the polarization normal to the crack, but weakly with the polarization parallel to the crack (Figure 8d).
We apply the Hilbert-Huang Transform (HHT) to analyze temporal variations in the CCCFs. This technique decomposes the time series data into several intrinsic mode functions (IMFs) through iterative Hilbert Transform (Huang & Wu, 2008). This technique can be regarded as an alternative to the time series filtering, but without concerns of phase shifting. We resample the temporal anisotropy data to one point per day and apply the HHT to analyze the response in different IMFs. The Piecewise Cubic Hermite Interpolating Polynomial is used to deal with glitch and irregularity of the input data. There are generally about 6-7 IMFs produced from the original anisotropy time series data, representing temporal change of seismic anisotropy in different time scales. We focus on the medium-term (seasonal) and the long-term (years) changes hereafter due to limited temporal resolution (∼4-month) as a result of 40-event stacking.

Medium-Term Variations
While the temporal change of seismic anisotropy seems to be affected by different mechanisms, we find that there is a clear seasonal variation in A ani and T 0 (Figures 9a and 9b). To understand the significance of seasonal response, we calculate the energy contribution of the seasonally varying IMF to the total energy. The calculation of energy contribution is defined as ∫ IMF 2 dt∕E , where E is the summation of every ∫ IMF 2 dt . We find the seasonally-varying IMFs in Figures 9a and 9b contribute the highest energy fraction of about 29%-56% ( Figure S6; Table S1 and S2). We find the seasonal variation in the anisotropy parameters is correlated with precipitation ( Figure 9c). HHT from A ani and T 0 show substantial fluctuations with amplitudes of about 2% and 0.001 s inferred from the seasonally-varying IMF, respectively (Figures 9a and 9b). The fluctuations reported here are robust as they compare A ani and T 0 at different time points, but only above the threshold of A ani measurement (i.e., greater than 6.25%). The periodicity of these fluctuations ranges from about 9 to 12 months (Figure 9d). We do not include the IMF analysis on the FSP time series since it does not display such a strong seasonal response. Figure 10 shows the statistics and cross-correlation between A ani and precipitation. In particular, we find that 86% of the precipitation between 2007 and 2010 are within the encouraging phases, during which the increase of seismic anisotropy is faster than the background trend. Only 14% of the precipitation events occur during the discouraging phases (Figure 10a), where the decrease of seismic anisotropy is faster than the background trend. Similar trends can be observed between T 0 and precipitation (Figure 10c), indicating a strong correlation between the isotropic velocity (i.e., 1/T 0 ), the A ani , and rainfall events. Overall, the seasonal variations coincide with the precipitation data where strong rainfalls during the wet season are closely associated with the increase of A ani and T 0 . These two parameters are also observed to decrease during the drought season.
Although the temporal resolution in our study does not permit a robust detection of temporal changes at shorter periodicity, seasonal change in 1/T 0 identified here corroborates the result from ambient-noise cross-correlation on the data at TCDPBHS (Hillers et al., 2014). Based on the observation of seasonal variation in anisotropy correlating with precipitation, we analyze the time lag between anisotropy change and the precipitation events. The cross-correlation between the anisotropy parameters and the precipitation is carried out through the calculation ∫ 0 ( ) ( − ) in which a(t) and p(t) are the anisotropy, (Figures 9a and 9b) and the precipitation time series (Figure 9c), respectively; T is the overall time span of 2007-2010 and is the time lag. The cross-correlation result is shown in Figures 10b and 10d, indicating there is a positive time lag (i.e., the anisotropy response lags the precipitation) of about 14 and 20 days for T 0 and the A ani , respectively. To test the robustness of these time lags, we repeat the measurements with 30-CCCF stacking and reconstruct the time series. While there is a similar time lag between A ani and precipitation of about 20 days, the seasonal signal is less apparent in T 0 , making it difficult to estimate the time lag.
This time lag might be indicative of the hydrological response of the fault zone to surface hydrologic events. We use the 20-day time lag to estimate the apparent hydraulic diffusivity D a ≈ h 2 /4t, where h is the depth, and t is the time lag between A ani and the precipitation. Assuming the depth h of 1,000 m, we obtain D a of about 0.14 m 2 /s. Uncertainty estimation from anisotropy curves derived from different numbers of CCCF stacking suggests the time lag is in the range of tens of days, resulting in D a on the order of 10 −1 m 2 /s. In general, the hydraulic diffusivity of rocks falls between 1 and 10 −4 m 2 /s (H. F. Wang, 2000). While in-situ hydraulic experiment in the Chelungpu fault suggests a very low diffusivity of about 10 −5 m 2 /s in the fault core area (Doan et al., 2006) where materials exhibit fine-grained, impermeable behavior , we conceive that the estimated hydraulic diffusivity here might be related to the fluid diffusion in the hanging wall side of the fault damage zone where microcracks are pervasive in the formations (H. Y. Wu et al., 2007). In particular, the fluid diffusion is more favorable in the Cholan fm., which consists of high porosity and permeable sandstones (T. -M. N. Chen et al., 2009). Our observations and the simple calculation above provide a first-order understanding of the hydraulic parameters in the fault damage zone. They may also potentially offer the opportunity to evaluate the competition between cohesion strengthening due to hydrothermal reaction and fault weakening resulting from increasing pore pressure (Tenthorey & Cox, 2006).
New observations of seasonal variations in seismic anisotropy presented here potentially help illuminate the mechanics responsible for the seasonal signal inside the fault zone. One possible mechanism is that the seasonal change in anisotropy parameters reflects the change of crack geometry or pore space in the fault zone. Shapiro and Kaselow (2005) showed that the change of differential stress may result in the change of pore space geometry, modifying the elastic moduli and seismic velocities. Pore pressure has been shown to play a significant role in extensional fracturing in low Δσ environments (Cox et al., 2010). It is conceivable that the Chelungpu fault-zone has low Δσ that favors cracks opening during precipitation events, which would produce changes in the anisotropy parameters such as A ani and T 0 . Such a small Δσ may be reconciled with the change of stress regime from thrust to normal-fault or strike-slip faulting after the rupturing of the Chi-Chi earthquake (Hung et al., 2009;W. Lin et al., 2007;Yeh et al., 2007).

Long-Term Fault-Zone Evolution
Beyond seasonal variations of the fault-zone anisotropy, a long-term trend can be obtained from the residual part of the HHT. The residual is obtained after all the oscillating IMFs are removed from the original time series. We found that the FSP direction shows a gradual rotation from 105° to 90° between 2007 and 2010, and stays  (Figure 11a). T 0 displays a slight increase over time, whereas the strength of anisotropy A ani shows a gradual reduction from 17% to 10%, and perhaps a slight increase after 2011 (Figures 11b  and 11c). We repeated the HHT analysis against the result from 10-event CCCF stacking and found the results are generally consistent ( Figure S7). There are slight differences in the anisotropy strength between 2011 and 2013, possibly due to differences in the strong short-term fluctuations.
Such a long-term trend in the anisotropy parameters is potentially indicative of fault zone evolution after the Chi-Chi earthquake between 2007 and 2013. We hypothesize that these observations reflect the continuation of the fault healing process over 10 years after the earthquake. Figure 11d illustrates a schematic model to discuss the healing process at different stages. During stage I (about 8 years after the 1999 Chi-Chi earthquake), vertical microcracks are likely oriented along 105° with anisotropy strength of about 17%. During stage II, the A ani reduces from 17% to 10%, where some microcracks are preferentially closed or sealed, with remaining cracks orienting at 90°. At stage III, the FSP direction remains relatively constant from 2011 to 2013, over 10 years after the earthquake.
It has been shown that the average S Hmax direction inferred from the borehole breakout in the Formation Micro Scanner (FMS) and the Formation Micro Image (FMI) analysis is about N115°E (H.Y. Wu et al., 2007). However, borehole breakouts in the fault damage zone of 950-1,100 m, just directly above the identified slip plane, deviate substantially toward fault-strike parallel (H.Y. Wu et al., 2007). Interestingly, W.  and Lin et al. (2010) also identified a similar change in borehole breakout directly above the slip plane in a neighboring hole (i.e., hole B), admittedly within a more restricted zone of ∼20 m. Perhaps, the S Hmax in the fault zone directly above the slip plane was fault-parallel 6 years after the Chi-Chi earthquakes in 2005, but changes rapidly to fault-perpendicular in 2007 and remains relatively stable since. Such a 90-degree switch in the S Hmax over time is not inconsistent with the change in paleostress over multiple seismic cycles (Hashimoto et al., 2015). While this is obviously a critical issue, at this stage, there is still inconsistency in the reported borehole breakout direction (e.g., Hung et al., 2009). Such a discrepancy remains to be resolved and we will refrain from further speculation.
The results of this study are subject to three possible sources of uncertainty: (a) borehole sensor rotation, (b) errors in earthquake location, and (c) variability in earthquake source depth. TCDPBHS uses the Galperin system to record ground motions. Lin et al. (2012) verified the sensor correction has accuracy with azimuthal uncertainty of about 5°. We minimize errors due to borehole sensor rotation by using a rotation increment of 5°. The error from earthquake locations would only contribute to our interpretation of the SWS results, since it indicates the overall anisotropy behavior along the ray path. We confine the earthquake incidence with a reasonable range of 40°, which accommodates location uncertainty. Finally, to accommodate variability in earthquake source depth, we use earthquakes in a narrow depth range between 9 and 14 km. This narrow depth range minimizes the effects of time delay in calculating anisotropy parameters. For CCCF analysis, we also constrain earthquakes to a confined incidence because previous workers have pointed out coda energy flux may vary with respect to the earthquake source location leading to an unstable CCCF (i.e., Emoto et al., 2015). Our selection of earthquakes minimizes any instability in CCCF analysis.

Conclusions
We present in-situ monitoring of seismic anisotropy using the 1-km deep, TCDP borehole seismic array near the Chelungpu fault, the main rupture fault of the 1999 Chi-Chi earthquake. We used shear-wave splitting (SWS) and coda cross-correlation (CCCF) to study seismic anisotropy in the shallow crustal region and in the fault damagezone. Although the SWS result suggests relatively weak anisotropy in the shallow crust, strong anisotropy in the fault-zone is detected through CCCF analysis. The fault-zone anisotropy revealed from CCCF analysis shows an average FSP direction of about 93° and a strong anisotropy of about 12% during the observational period between 2007 and 2013. Although the fault-zone anisotropy could result from both stress-induced and intrinsic anisotropy, pervasive vertical microcracks after the M7.6 Chi-Chi earthquake likely dictate the observed anisotropy parameters. The fault-zone anisotropy also displays a seasonal periodic pattern of about 9-12 months and is correlated with rainfall events. These seasonal variations indicate the deep fault (∼1 km) is sensitive to surficial hydrologic processes. The multi-year evolution of fault-zone anisotropy shows three key observations: (a) a gradual rotation of the FSP direction by about 15° during the first 4 years of observation, (b) a notable reduction of anisotropy strength from 17% to 10%, and (c) an increase in the isotropic shear-wave velocity. These combined observations indicate slow, but continuous healing of the Chelungpu fault due to crack closure. This study highlights a rare opportunity to interrogate the temporal evolution of the fault damage zone stress state in-situ, and to illuminate the response of the damage zone from external modulation (i.e., precipitation) during fault healing.

Data Availability Statement
The TCDP seismic data and catalog may be found at https://data.mendeley.com/datasets/pxmt37cr9p/1. Details are described in Supporting information. The precipitation data is available from the Central Weather Bureau, Taiwan (https://e-service.cwb.gov.tw/HistoryDataQuery/). The software MSNoise is an open source code at http://www.msnoise.org/(last access July 2021). The HHT was calculated via MATLAB 2019a.