Monitoring In-Situ Seismic Response on Rock Slopes Using Ambient Noise Interferometry: Application to the 2019 Changning (Mw 5.7) Earthquake, China

Study of the mechanical response of rock slopes to moderate earthquakes is important for understanding the local rheology of landslide and earthquake interactions and for mitigating the risks associated with subsurface geological processes in tectonically active mountainous belts. To complement existing point measurements from surface observations (e.g., global positioning system and interferometric synthetic-aperture radar measurements), measuring the ambient noise-based velocity change ( δ v / v ) allows for remote observations of mechanical state changes of the slope, at depth and continuously in time. We herein investigate the seismic responses of the Pubugou rock slope, a typical steep rock slope in south-west China, to the 2019 Mw 5.7 Changning earthquake. We apply ambient noise interferometry to the slope and measure the coda wave velocity changes at frequencies from 2 to 20 Hz with a 1-h temporal resolution, 2 days before and 14 days after the earthquake. We observe a significant co-seismic wave velocity decrease caused by the Changning earthquake of up to 0.9% followed by a gradual logarithmic recovery process over 2 weeks. The earthquake-induced stress sensitivity of δ v / v on the slope is estimated as ∼ 3.2 × 10 − 8 Pa − 1 . Through the analysis of the co-seismic and post-seismic δ v / v with different time lapses of the coda, we characterize the healing process on the slope and also constrain such changes to 75 m in depth. This study highlights the possibility of quantitatively characterizing the slope weakness using moderate earthquakes in mountainous areas in the future.


INTRODUCTION
Deep-seated rock slopes are widely distributed in the mountainous areas of Sichuan Province, which is one of the most active geohazard regions in China. There are over 300,000 slopes in this area, which are susceptible to a high level of landslide activity (Lin and Wang, 2018). Simultaneously, it is also a seismically active area. Increased seismicity has been recognized in the last decade since the 2008 Wenchuan earthquake (Chigira et al., 2010). According to the earthquake catalogs from Sichuan Earthquake Administration, over 4,000 earthquakes occurred in 2019 with magnitudes ranging from Mw 1 to Mw 5.7.
Mass movement on the damaged rock slopes constitutes a major geological hazard, damaging infrastructure such as dams, roads, railways, and bridges and leading to loss of life. Among a variety of physical parameters (e.g., atmospheric pressure, tide, temperature, and rainfall), seismicity has a strong impact on the damage evolution of the rock slopes as dynamic strain is applied by seismic shaking. Previous studies on seismic hazard assessment on slopes after earthquakes have been based on the statistical analysis of regional inventories of earthquake-triggered landslides (Keefer, 1984;Keefer, 2002). Due to a lack of in-situ data and a cost-effective method that is sensitive to the changes of elastic properties, there have been few quantitative studies on earthquake-induced temporal changes on slopes.
Observational evidence (Larose et al., 2015) and theoretical models (Colombero et al., 2017) suggest the earthquake-induced damage is associated with changes in the material's elastic moduli, which lead to failure of the rock slopes. Owing to the heterogeneous elastic nonlinearity, it is hard to quantify such damage on the rock slopes at scales ranging from macroscopic fractures to microscopic contact changes between grains. Nevertheless, it is widely recognized that loss of rigidity is a fundamental signature of the damage development. Therefore, this makes monitoring the changes in elastic wave velocity (Murnaghan, 1951) an ideal method to remotely assess the internal damage development of the rock slopes.
Among the emerging techniques in this field, one promising monitoring approach is ambient seismic noise cross-correlation. This is a passive technique that enables retrieval of impulse responses through cross-correlation of ambient seismic noise recorded at any two sensors. Depending on the multiple scattered coda waves retrieved by the ambient crosscorrelation technique (Larose et al., 2006), it is possible to further monitor stress changes of the medium by inferring velocity changes from the phase shift in the coda at different times (Snieder, 2006). This forms the basis of ambient seismic noise interferometry for time-lapse applications. Ambient seismic noise interferometry has been used in geophysics for more than 30 years, e.g., to study the dynamic evolution in faults by monitoring velocity changes caused by nearby earthquakes (Brenguier et al., 2008). This method has been proposed to study the precursor instabilities responsible for landslides. A significant velocity reduction was reported several days before the failure on the Pont Bourquin clay soil of the Swiss Alps (Mainsant et al., 2012). Moreover, the method enhances the knowledge of deformation by environmental variations (e.g., seasonal fluctuations and groundwater infiltration) on slopes with materials composed of clay soil and/or volcanic deposits (Larose et al., 2015).
Monitoring changes of rock slopes during earthquake shaking and the subsequent recovery phase is used to study the dynamic elastic properties of the rock slopes and also to evaluate the damage level. However, to the best of the authors' knowledge, only few observations with limited typology of the slope have successfully tracked co-seismic changes. Even less understood is the healing process due to post-seismic changes when moderate dynamic stress solicitation is applied on rock slopes by a distant earthquake. This paper focuses on the co-seismic and post-seismic changes on the Pubugou rock slope due to the Changning Mw 5.7 earthquake in 2019. The coda wave velocity changes with a 1-h temporal resolution are measured by applying ambient seismic noise interferometry to continuous recordings on the slope. With these results, some possible explanations and implications of the observations are discussed.

Study Area
As shown in Figure 1A, the Pubugou rock slope is in the middle of the deep-seated bare bedrock alpine valleys of Dadu River between the west margin of Sichuan Basin and the Tibet Plateau. The rock slope is ∼800 m from the dam of the Pubugou Hydropower Station (3,600 MW capacity) and is surrounded by active faults characterized by complex Cenozoic structures. There are intense deformations and high levels of seismic activity in this area. Figure 1B illustrates that the east-facing slope with its height of about 1,180 m above mean sea level (m a.m.s.l.) consists of highly weathered tuffaceous rocks and diabases, in addition to gravels with thin clays on a shallow layer at the upper rear part. The gradient of the slope is approximately 45°. The surface displacement monitoring indicates progressive toppling of the upper part (980-1,180 m a.s.m.l.) at an annual displacement rate of 21 mm/yr. The damaged lower part (< 980 m a.s.m.l.) has been reinforced by removing shallow weathered gravels and rocks and inserting ∼100 anti-slide piles with cement grout into drilled holes. Thus, the annual displacement rate in this area is now approximately 3 mm/yr. According to the drilling data provided by a geology study (Yang et al., 2013), intact rocks exist at an average depth of ∼70 m. The drilling samples also suggest that the weathering level on the rocky material varies with depth. The strongest weathered material appears at a shallow depth of a few meters.

Data
As illustrated in Figure 2A, two observational huts (T01 and T02) were located ∼157 m apart on either side of the center part (∼990 m a.s.m.l.) of the slope. In such way, a global slope coverage due to limited sensors as well as strong mobile network signal condition was expected. Each hut had a threecomponent short period seismometer (GL-CS2, GEOLIGHT ™ ) installed, which allowed a relatively flat frequency response from 0.5 to 50 Hz. Each seismometer was digitized by a 24bit EDAS-24 GN (GEOLIGHT ™ ) at a 100 Hz sampling rate with a system clock error less than 10 − 6 s. The seismic noise data were stored locally in 1-h long records and simultaneously transmitted to the server at Sichuan Earthquake Administration via a mobile network.
A bi-frequency global positioning system (GPS) receptor (GMX902, Leica ™ ) was installed on the roof of each hut, together with a GPS antenna at each hut, to measure the surface deformation along the central part of the slope since April 24, 2019. The hourly global navigation system satellite Frontiers in Earth Science | www.frontiersin.org January 2021 | Volume 8 | Article 610181 (GNSS) measurements were transformed into daily values by averaging over every 24 h. Daily rainfall data were gathered from a weather station on the east bank of the rock slope. On June 17, 2019, an Mw 5.7 earthquake occurred in Changning county, a region of Sichuan Province in southwest China. Despite the distance between the rock slope and the epicenter of the earthquake being ∼221 km, the seismometers on the slope clearly recorded earthquakeinduced ground motion. In Figure 2B, the seismic shaking intensity is estimated by directly computing the peak ground velocity (PGV) for its vertical component. The measured PGV is 2.34 cm/s for T01 and 2.96 cm/s for T02. Hence, the

METHOD
The hourly ambient seismic noise was pre-processed at each seismometer following a standard routine procedure (Larose et al., 2015): the hourly seismic noise was first normalized in the frequency domain between 2 and 20 Hz using the whitening method and then normalized in the time domain using the clipping method with amplitudes exceeding three standard deviations. The purpose of the above steps was to enhance the specific frequency band of the ambient noise and also reduce the impact of spatially isolated noise sources such as earthquakes. Finally, the vertical components of the pre-processed noise values between the seismometers were cross-correlated for each hour with a 15 s time lag. As the cross-correlation functions were asymmetric, due to the anisotropy of the noise propagation direction, the causal and acausal parts of the cross-correlation functions were further averaged. Figure 3B illustrates the 384 hourly correlograms h i (i denotes the hour during the period of interest from June 15 to June 30) calculated, along with the reference h ref ( Figure 3A) that was obtained by averaging all the correlograms.
To understand how changes in seismic velocities can be measured from phase shifts in the correlograms, consider a ray that travels a distance L between two seismometers in time t at velocity v. These quantities are related via

L vt
(1) Taking the differential of L, For a homogeneous velocity change, the ray path does not change (δL 0). This yields a relation between the velocity change δv and the change in arrival time δt: The change in arrival time scales linearly with lag time, resulting in stretching or compression of the waveform following a decrease or an increase, respectively. Thus, to quantify the temporal changes, the stretching method was applied, which consisted of a grid search testing over a series of candidate velocity changes δv/v k between h i (t) and the stretched reference signal h ref (t(1 + δv/v k )) within a given time window [t1, t2]: The apparent velocity changes δv/v can be obtained by maximizing the correlation coefficient CC(δv/v k ). The uncertainty on δv/v, which is the indicator of the measurement reliability, can be estimated using the theoretical formula proposed by Weaver et al. (2011): where T is the inverse of the frequency bandwidth and ω c is the central frequency. The lower the uncertainty is, the robuster the measurement is.

Co-Seismic Change
A time window from 1 to 12 s was used ( Figure 3) to measure the in-situ global temporal change using scattered coda waves. Figure 4A illustrates the general evolution of the apparent velocity change together with its color-coded uncertainty as a function of time during the period of interest from June 15 to June 30. The uncertainty (∼ 10 − 5 ) is one or two orders of magnitude smaller than the apparent velocity change (∼ 10 − 3 ), suggesting that the measurement is consistently good in quality.
We observe a rapid co-seismic velocity drop by ∼0.9% on the slope. Simultaneously, the recorded PGV is 2.34 cm/s for T01 and 2.96 cm/s for T02 due to the Changning earthquake. Thanks to the measured PGV, the earthquake-produced dynamics stress on the slope can be estimated using the following equation (Hill et al., 1993): where G is the bulk modulus and v p is the bulk velocity. As with the estimation of physical properties of tuffaceous rocks (Yang et al., 2013), we use v p 1.5 km/s and G 1.6 GPa. Hence, Δσ d is approximately ∼280 kPa, and a velocity-stress sensitivity of ∼ 3.2 × 10 − 8 Pa − 1 can be estimated. The seismologists (Peng et al., 2010;Gonzalez-Huizar et al., 2012) have reported that even though the earthquake generating PGV remains as low as millimeters per second at a distance of thousand kilometers, the induced stress loading still exists at least 10,000 Pa and has the potential to remotely trigger seismicity. Therefore, such a strong but sudden decease in velocity cannot be influenced by the undrained loading effect due to the groundwater-induced liquefaction of the slope (Lecocq et al., 2017), or by the thermal loading effect due to the air temperate change (Tsai, 2011). The only possible mechanism is attributed to the opening of new or pre-existing cracks/fractures on the slope which decrease the elastic modulus due to the stress loading of earthquake shaking (Bontemps et al., 2020). We interpret that the observed co-seismic velocity drop is caused by the weakness of the slope due to the earthquake-induced dynamic stress transient.

Post-Seismic Change
Almost immediately after the co-seismic velocity drop, a clear logarithmic recovery of the velocity change to the pre-earthquake state is observed over approximately 14 days. Such a healing increase back toward its initial value following a logarithmic evolution after moderate solicitation is referred to as "slow dynamics" (TenCate, 2011).
Slow dynamics is a nonlinear elastic response of the material after imposing a strain of moderate amplitude that does not generate any macroscopic damage and has been found to be universal in granular solids of various compositions ranging from that of the Earth's crust to the inter-grained microscopic cracks of sedimentary rocks or concrete samples of the order of micrometers (TenCate et al., 2000). Hence, we interpret this recovery phase as a nonlinear response to the re-compaction of the opened fractures and micro-cracks due to seismic shaking. Note there was an earthquake sequence that began with an Mw 5.7 earthquake on June 17 and comprised ∼600 aftershocks for more than 2 weeks. Within this context, several moderate earthquake-induced ground motion events were recorded on the slope during this period. Among them, there was a 0.1% co-seismic velocity decrease caused by the largest Mw 5 aftershock on June 22 during the healing process, with a PGV value of 1 cm/s on the slope.
In addition, a small velocity fluctuation was observed at the end of the recovery phase. It was noted that such velocity fluctuations occurred after 2 days of precipitations with a total of ∼90 mm rainfall. Rather than earthquake shaking (with loading effects that have low PGV values (< 0.1 cm/s)) during these days, the possible main mechanism response for such velocity fluctuations was rigidity decrease due to an augmentation of the water content/pore saturation in the highly weathered granular material of the slope above a certain rainfall threshold (Guzzetti et al., 2008). However, the complex interactions between seismic shaking and rainfall are out of the scope of this work and will be demonstrated in a separate paper.
To quantify the healing process, the logarithmic evolution of δv/v can be fitted as follows: δv v A − m log 10 (d) where d is the recovery time in hours, A is the extrapolated δv/v after the co-seismic shaking, m is the slope of the logarithmic decay, and 10 A/m represents the recovery time. Figure 4B shows the estimation of these parameters with (June 18-27) and without (June 18-22) the Mw 5 aftershock on June 22. The characterized recovery times (10 A/m ) are 146 and 177 h, respectively.
This suggests that the rigidity of the slope experiences a faster recovery process with the effect of additional small seismic event(s). It also implies that even seismic shaking with its PGV value as low as 1 cm/s can accelerate the healing recovery processes by altering the slope's weakness. A possible mechanism for this effect is the re-arranging of existing grains/ cracks that favor the closing of opened earthquake-generated micro/macro-fractures (Bontemps et al., 2020).
It is worth noting that compared with ambient noise interferometry, the horizontal displacement ( Figure 4C) measured at each seismometer site did not show significant co-seismic and post-seismic effects. One reason is that due to low temporal resolution (1 day) of the GNSS measurements, it was hard to track sudden earthquakeinduced changes over such a short time. However, such measurements hardly reveal small changes of mechanical properties of the material because they are less sensitive to the state of stress, rigidity, or damage due to their surface measurement configuration.

Different Times in the Coda
The changes are assumed to be in a spatially global homogeneous medium (Snieder, 2006). Thus, velocity change is proportional to a time shift in the later arriving coda waves. However, as a strong heterogeneous material, velocity change in the rock slope is not global, but is influenced by the time in the coda with respect to the sampled domain. Therefore, instead of analyzing a long coda segment as in the previous section, the stretching technique is applied to a series of sliding windows of 3 s in the coda, with the window width corresponding to six periods of applied lowest frequency (2 Hz). Figure 5 illustrates the velocity changes within seven consecutive shorter time windows centered around t c 1.5, 3, 4.5, 6, 7.5, 9, and 10.5 s. The strongest co-seismic velocity reduction and post-seismic recovery process were found at the earliest time of t c 1.5 s in the coda. With increasing time, the magnitude of the co-seismic velocity drop decreased from ∼0.9% to ∼0.1%. Simultaneously, a weaker logarithmic post-seismic recovery process in the later arriving coda was observed. The weakest coseismic and post-seismic velocity perturbations were observed during the last time window of t c 10.5 s in the coda.
As the early coda waves are mostly sensitive to the changes at shallow depths, as well as the weaker bonds between the elements of the weathered tuff rocks, we suggest that earthquake-induced strong ground motion can introduce more opening cracks/ fractures, which is the signature of mechanical damage in the form of crack appearance, at a shallow depth of the higherweathered rock material in the slope.
It is worth noting that the sensitivity to absolute depth requires detailed measurements of the scattering properties of the slope, which are not available for the study area. Nevertheless, the relative depth resolution of coda waves measured at different time lapses of the coda using existing theoretical and numerical testing can be discussed.

DISCUSSION
The in-situ apparent velocity evolution on the rock slope due to moderate earthquake shaking raises at least two questions. First, what is the depth sensitivity of the measurements? Second, why is the slow dynamics effect due to earthquake shaking related to the damage features?
To answer the first question, we estimate that the first-order depth resolution for the velocity decrease is due to earthquake. Although the depth sensitivity of coda waves has not yet been fully solved (Obermann and Hillers, 2019), as the study area is neither a layered medium nor a completely solid rock material, we thus evaluate the bulk sensitivity of the scattered coda waves by considering a two-dimensional body wave sensitivity kernel formulation (Obermann and Hillers, 2019): where s and r are the positions of the stations, r 0 is the local velocity variation, and t is the lapse time. The energy propagator p is calculated by the radiative transfer solution approximation for isotropic scattering of the medium: where H(x) is the Heaviside function and c is the energy velocity.
With an empirical scattering mean free path of ∼10 m as well as a bulk energy velocity of ∼1,100 m/s, Figure 6A gives the theoretical depth sensitivities with respect to the different time lapses (t c ) of the coda, in which sensitivity is normalized by a body wave kernel. Figure 6B plots the normalized depth sensitivities with values greater than 10% as a function of measured co-seismic velocity reduction at each lapse time. Figure 6B shows that the δv/v measurement is most sensitive to changes within the top ∼100 m, and velocity reduction varies with its depth above ∼75 m. This result is in agreement with geological investigation that the intact For the second question, due to their direct sensitivity to the higher-order elastic properties of materials, scattered coda waves are more sensitive to the nonlinear elastic response of the material than direct waves (Xie et al., 2018). Compared with nondestructive testing techniques based on common elastic waves, the emerging methods based on coda waves have shown great promise in terms of being more accurate and sensitive in detecting the initiation of damage in various solid materials as early as possible (Xie et al., 2019). As a particularly attractive nonlinear elastic effect, time-logarithmic recovery (slow dynamics) back to the unperturbed elastic modulus following a sharp drop in the elastic modulus (fast dynamics) has been reported in response to moderate mechanical solicitation. Such nonlinear responses were first probed with nonlinear resonant ultrasound spectroscopy and can now be monitored efficiently via velocity changes from coda wave interferometry at various scales (Tremblay et al., 2010).
As far as earthquake-triggered landslides are concerned, some studies have suggested that in contrast to the postseismic recovery during the dry season, rainfall-produced fluid activities disturbed the post-seismic relaxation process causing a slower recovery time (Bontemps et al., 2020). Laboratory studies have also suggested that the characterized  Frontiers in Earth Science | www.frontiersin.org January 2021 | Volume 8 | Article 610181 7 recovery times are dependent on the damage of the medium (Tremblay et al., 2010).
In addition, in contrast to seismometer's record, daily GNSS has been applied comprehensively to measure subtle surface mass deformation at a high confidence level of millimeter scale with 24-h static data. It is presumed that the deformation continuously accumulates from months to years, where the GNSS suffices in such high precision level. Nevertheless, there are indeed strong-motion scenarios with rapid deformation (e.g., storm surge loading, pre-eruption volcanic unrest, and earthquakes), which require epoch-wise GPS displacements to capture motions on a wide sub-daily timescale from seconds to hours. However, high-rate GPS solutions to detect sub-daily deformations in a reliable manner are still of great challenges as the magnitudes of such signals are usually close to the lower bound of GNSS carrier-phase measurement precision (Bilich et al. (2008); Geng et al. (2017)). Therefore, so far, it is hard for current GNSS techniques to reveal such small changes of mechanical properties of the rock slope due to their low sensitivity.
To summarize, from scientific point of view, this work quantitatively characterizes the weak but solid dynamic loading effect due to a distant earthquake and subsequent time-dependent recovery process in terms of elastic change using ambient noise interferometry. It is the first time to reveal such physical process inside the rock slope at such highly temporal (1-h) and stresssensitive (∼ 3 × 10 − 8 Pa − 1 ) ways due to seismic shaking at a distance over 200 km. It facilitates the understanding of on-site damage evolution of the rock slope during the earthquake. From engineering point of view, since seismic shaking is ubiquitous in tectonically active mountainous belts, it has the potential to quantitatively characterize the slope weakness, which is susceptible to nonlinear elastic changes in terms of velocity changes by monitoring the post-seismic relaxation process. It facilitates the in-situ seismic hazard assessment of the rock slope during the earthquake which was previously based on the statistical analysis of regional inventories of earthquake-triggered landslides.

CONCLUSIONS
This study has applied in-situ ambient noise interferometry to two seismic stations installed on a slope and has measured coda wave velocity changes at frequencies between 2 and 20 Hz with a 1-h temporal resolution, 2 days before and 14 days after an earthquake. The findings are as follows: (1) Co-seismic wave velocity decreases caused by the Changning earthquake of up to ∼0.9% were followed by a gradual logarithmic recovery process over 2 weeks. An earthquakeinduced stress sensitivity of δv/v on the slope was estimated as ∼ 3.2 × 10 − 8 Pa − 1 .
(2) By analyzing the co-seismic and post-seismic δv/v with different coda time lapses, the healing process on the slope was characterized, and such changes were constrained to ∼ 75 m in depth.
(3) It may be possible to quantitatively characterize slope weakness using moderate earthquakes in mountainous areas in the future.

AUTHOR CONTRIBUTIONS
HH and FX designed the experiments, analyzed the results, and wrote the paper. SD collected the earthquake catalogue and processed the data acquired. All authors contributed to the manuscript revision.