S Coda and Rayleigh Waves From a Decade of Repeating Earthquakes Reveal Discordant Temporal Velocity Changes Since the 2004 Sumatra Earthquake

Temporal changes in the subsurface seismic velocity structure reflect the physical processes that modulate the properties of the media through which seismic waves propagate. These processes, such as healing of the surface damage zone and deep crustal deformation, are described by similar functional forms and operate on similar timescales, making it difficult to determine which process drives the observed changes. We examine earthquake‐induced velocity changes using the measured lag‐time time series τ(t) of the repeating earthquake sequences since the 2004 Mw 9.2 Sumatra and 2005 Mw 8.6 Nias earthquakes. The S coda velocity changes (δVS, equivalent to −τS) recover steadily during the 2005–2015 period. The Rayleigh wave velocity changes (δVLR, or −τLR) undergo transient recovery, followed by a strong δVLR reduction in late 2007. δVS recovery is most likely driven by deep processes, whereas the temporal breaks in δVLR recovery in 2007 mostly reflect surface damage and healing induced by the strong ground motions of the 2004 Mw 9.2, 2005 Mw 8.6, 2007 Mw 8.4, and Mw 7.9 Bengkulu and 2008 Mw 7.3 Simeulue earthquakes. The observed differences between the temporal variations in δVS and δVLR can distinguish deep processes from healing of the surface damage zone.


Introduction
Long-term monitoring of temporal seismic-velocity variations in the crust is an important and long-sought goal in geophysics because these temporal changes can be proxy measures for the mechanical processes and timescales of crustal responses to earthquake slip. Numerous studies have suggested that these temporal velocity changes are often associated with the peak ground velocity (PGV) or acceleration (PGA) that is induced by the strong ground motion (GM) of an earthquake (Hobiger et al., 2016;Rubinstein & Beroza, 2004a;Rubinstein et al., 2007;Schaff & Beroza, 2004;Takagi et al., 2012;Wegler et al., 2009;T. Yu & Hung, 2012;Y.G. Li et al., 2006). Strong GM damages the subsurface sediment layer, which is subjected to low confining pressures of only a few tens of MPa (L. Li et al., 2017;Sawazaki et al., 2015;Wu et al., 2009;Zhao & Peng, 2009). The subsequent velocity recovery, which probably reflects healing of the subsurface damage zone, is typically characterized by a logarithmic time dependence (so-called "slow dynamics") (Brenguier et al., 2008;Peng & Ben-Zion, 2006;Snieder et al., 2017;TenCate et al., 2000). Similar pattern of temporal velocity changes may also be associated with aseismic slip (or slip rate), postseismic afterslip (Rivet et al., 2011;Yu et al., 2013a) and viscoelastic relaxation in the bulk crust. Postseismic transients due to afterslip and dilatancy recovery occur adjacent to the earthquake slip zone and can modulate the crack density and pore pressure, yielding recoveries in S wave velocity and Poisson's ratio (Fialko, 2004;Fielding et al., 2009;Peltzer et al., 1996). Viscoelastic relaxation mainly occurs in the lower crust and upper mantle over longer timescales and can restore the stress perturbation induced by coseismic slip (Barbot & Fialko, 2010;Qiu et al., 2018;Tang et al., 2019). Postseismic transient processes due to dilatancy recovery (Fialko, 2004;Jónsson et al., 2003), afterslip (Froment et al., 2013), and viscoelastic relaxation all follow similar logarithmic decay patterns. However, it is difficult to unambiguously associate an observed logarithmic velocity recovery with a specific physical mechanism and/or depth (Obermann et al., 2016(Obermann et al., , 2013Yang et al., 2019) since all of the processes that influence the subsurface damage zone and postseismic transients operate on similar timescales. For example, both the healing of the fault zone (Y. G. Li et al., 1998) and decay of postseismic deformation (Fialko, 2004) are prominent over several years after the 1992 M w 7.3 Landers earthquake.
W. Yu et al. (2013bYu et al. ( , 2013a detected repeating earthquakes (REs) near the Sumatra Subduction Zone and used these repeaters to probe the temporal changes in seismic velocity induced by the 2004 M w 9.2 Sumatra and 2005 M w 8.6 Nias earthquakes during [2005][2006][2007][2008]. They found that the temporal velocity changes in the high-frequency (HF, 0.5-2.0 Hz) S wave codas (δV S ) and long-period (LP, 0.03-0.1 Hz) Rayleigh waves (δV LR ) followed similar logarithmic recovery timescales after the 2004 and 2005 earthquakes that were potentially linked to strong GM-induced surface damage or/and postseismic afterslip. However, the ability to distinguish the velocity changes induced by deep processes from healing of the subsurface damage zone remains a challenge, as mentioned above.
Cracks produced by a previous earthquake can multiply within the subsurface damage zone as a result of the strong GM induced by subsequent earthquakes (Rubinstein & Beroza, 2004a). If the logarithmic velocity recovery reflects healing of the subsurface damage zone, then the redamage process from subsequent large earthquakes could result in further reduction of the seismic velocity, reversing the established recovery pattern (Rubinstein & Beroza, 2004a;Sawazaki et al., 2018;Vidale & Li, 2003). However, velocity recovery is likely to continue without reversal if it is driven by postseismic afterslip or other deep processes, even in the event of another earthquake. Therefore, the ability to detect such a bifurcation in the temporal velocity recovery patterns would help discriminate between the healing of surface damage zones and slow deformation in the crust.
The 2004 M w 9.2 Sumatra earthquake (Ammon et al., 2005;Chlieh et al., 2007;Park et al., 2005;Subarya et al., 2006) has been followed by a series of large earthquakes that have occurred in the proximity of the 2004 rupture zone and in the Wharton Basin, northeast Indian Ocean (Figure 1; see Table S1 in the supporting information for the event IDs). Examples of the former include the 2005 M w 8.6 Nias earthquake (Hsu et al., 2006;Ishii et al., 2007;Konca et al., 2007Konca et al., ), 2007 and M w 7.9 Bengkulu double earthquakes (Konca et al., 2008;Tsang et al., 2016), 2008 M w 7.3 Simeulue earthquake (Morgan et al., 2017), and 2010 M w 7.8 Banyak earthquake, and examples of the latter include the 2012 M w 8.6 and M w 8.2 double earthquakes H. Yue et al., 2012;Wei et al., 2013). The strong GM generated by these large earthquakes can repeatedly damage the surface layer and interrupt healing, whereas the large postseismic afterslip from the 2004 M w 8.8 Sumatra and 2005 M w 8.2 Nias events (Chlieh et al., 2007;Feng et al., 2015;Hoechner et al., 2011;Hsu et al., 2006) (Figure 1), among other slow deformation processes, potentially drives continuous velocity changes. The main goal of this work is to examine how seismic-velocity change measurements from different attributes can help differentiate between these concurrent processes.
We extend our previous work (W. Yu et al., 2013a) by incorporating seven additional years of regional seismic waveform data, which expands our RE data set to the 2005-2015 time period. We discuss the data processing procedures and summarize the observed temporal velocity changes in δV S and δV LR in the following sections, where we highlight a clear bifurcation in the δV S and δV LR patterns: δV S follows a monotonic logarithmic recovery, while δV LR deviates from the background recovery rate in late 2007. Through a forward modeling approach, we explore a series of perturbed velocity models and examine the sensitivity of these models against the observed attributes, to reconcile the first-order pattern shown in the observed lag-time time series τ(t) of the HF S coda and LP Rayleigh waves.  (Table S1). The repeating earthquake (RE) sequences are represented by small stars, and the background seismicity is represented by the gray dots, all of which is superimposed on the modeled coseismic slip (contours) and postseismic afterslip (colored area) of the 2004 Sumatra (Chlieh et al., 2007Nias (Hsu et al., 2006 (Table S2) are indicated by small yellow stars, whereas the REs with high cross-correlation (CC) coefficients and no further relocation analysis are indicated by small white stars; the source parameters of the latter are not provided. The REs are located near the afterslip zones of the 2004 and 2005 events (Chlieh et al., 2007). The solid and dotted black lines indicate the trench (Bird, 2003) and slab depths at 50-km intervals (Gudmundsson & Sambridge, 1998), respectively. The seismic and GPS stations used in this analysis are indicated by the open triangle and open squares, respectively. An enlarged map near the Nias, Banyak, and Simeulue Islands is provided to better present the 2005 Nias earthquake RE locations, with the dotted gray box indicating the location of the displayed region in the main figure.

REs Recorded by Seismographic Station PSI
The procedural details for the RE detection and relative location assessment within a given earthquake sequence have been discussed in previous studies (Wen, 2006;W. Yu, 2013;W. Yu et al., 2013b;W. Yu & Wen, 2012). Several previously identified RE sequences continued to occur in the afterslip zones of the 2004 and 2005 earthquakes until the end of 2015 (Figure 1). The high-precision relative hypocenter locations of these REs (relative to the estimated rupture area of the selected RE sequences) are presented in the supporting information (Text S1, Figure S1, and Table S2). Seismographic station PSI, which was installed in Parapat, Sumatra, Indonesia, in 1993, is equipped with a Streckeisen STS-2 broadband sensor. Station PSI detected the REs immediately after the 2004 and 2005 mainshocks, and it is in an advantageous position to potentially detect the δV induced by the 2005 event since the seismic waves from the RE sequences associated with the 2005 event propagated through the main slip and surface damage zones of the 2005 event as they traveled to station PSI.

τ(t) Measurements Using Coda Wave Interferometry
Coda wave interferometry is used to measure the τ(t) of the HF S coda and LP Rayleigh waves relative to the direct P wave onset in a given RE sequence (Poupinet et al., 1984;Snieder, 2006;Snieder et al., 2002). The HF waveforms are convolved using a zero-phase two-pole Butterworth band-pass filter with 0.5-and 2.0-Hz corner frequencies, and the LP Rayleigh waves are band-pass filtered between 0.03 and 0.1 Hz. The reference event of a given RE sequence is chosen as the last event before late 2007, when the onset of an additional increase in the Rayleigh wave lag times (τ LR ) is detected (to be discussed in section 3). We measure the differential time of the first 4-s HF P wave between the target and reference events via the cross correlation (CC) of each RE pair. We adopt 4-and 40-s sliding windows for the HF coda and LP surface waves, respectively, with a 95% overlap between successive time windows for both waves. The τ(t) for the HF S coda and LP Rayleigh waves are computed as the time differences between either the HF S coda or LP Rayleigh waves, and the HF P wave onset using a sliding-window CC approach. We emphasize that a 40-s time window is sufficient to include the entire cycle of a 20-s period Rayleigh wave packet. A 50-ms sampling interval of the original time series is interpolated to a 5-ms sampling interval to achieve subsample precision in time, which can alternatively be attained by interpolating the resulting CC time series. We find that the τ(t) calculated via CC of the interpolated time series is identical to that obtained via interpolation of the resulting CC time series. The τ(t) for the HF S coda waves (τ S ) can be expressed as follows: where t refers to the arrival times (or lapse times) of the initial 4-s HF P (subscripted "HF P") and S coda (subscripted "S coda") waves between the target (superscripted "trg") and reference (superscripted "ref") events of a given RE pair. The τ(t) for the 20-s dominant period long-period (LP) Rayleigh waves (τ LR ) can be expressed as follows: The τ(t) measurements can eliminate the uncertainties due to the relative origin time error of the event catalog and Global Positioning System (GPS) clock drift error between the two events in a given RE pair. τ(t) at the HF P wave onset (τ P ) is negligible due to media velocity changes because the waves are aligned relative to the direct P wave onset. Figures 2a and 2c Figures 2a and 2c, respectively. Signal-to-noise ratios (SNRs) are used to sort the time segments with sufficient signal strength and are defined as follows: The noise amplitude is measured as the maximum amplitude in the 5-to 15-s time window before the P wave arrival, and the signal amplitude is measured as the maximum amplitude in a given time window. W. Yu et al. (2013a) previously adopted SNR thresholds of 5 for HF S coda waves and 10 for LP Rayleigh waves. However, a minimum threshold of SNR = 10 with this new data set would exclude many τ LR measurements, yielding few temporal samples for the RE sequences. Therefore, we explore suitable SNR thresholds for the LP Rayleigh waves of smaller events. We find that a threshold requiring at least 80% of the 40-s lapse time segments to satisfy SNR of ≥5 and a CC coefficient of ≥0.9 between each waveform pair yields a reasonable number of τ LR measurements and stable τ LR values. We acquire 471 τ(t) measurements from the HF S coda and LP Rayleigh waves using the above-mentioned SNR and CC thresholds.

δV S and δV LR Measurements and Uncertainties
The observed τ S (t) for the 2005 Nias earthquake RE sequences is stretched in the later part of the coda waves (see section 3), similar to previous findings (Lobkis & Weaver, 2003;Poupinet et al., 1984). A fractional seismic-velocity change in the S wave coda (δV S ) of the RE pair can be approximated by a negative τ S (t) slope (dashed lines in Figure 2b).
where Δt S is the lapse time difference between the S wave coda at lapse time t and the S wave onset, with the corresponding lag time difference Δτ S (t). δV S is computed from the negative slope of time delays of the S codas that yields a path-averaged value (Poupinet et al., 1984;Rubinstein et al., 2007;Sawazaki et al., 2015;Schaff & Beroza, 2004). Linear regression is used to estimate the parameters of the slope b (−δV S ) and intercept a whose τ S values that meet the SNR and CC coefficient thresholds. The 2σ standard deviations in δV S are twice the computed root-mean-square (RMS) of τ i S −τ i S , where b τ i S ¼a þ bΔt S , which is the predicted value based on the inferred a and b values. However, τ LR is nearly The amplitudes of the filtered waveforms in each selected time window are normalized to the maximum amplitude in that window. The shaded regions correspond to the waveform segments with signal-to-noise ratios (SNR) of ≥5 and cross-correlation (CC) coefficients of ≥0.9, which are used to compute δV. The unstable and oscillatory τ(t) values at the beginning of each lower plot are due to low SNR values. The dashed black and light-red lines in (b) correspond to the −δV S values that are estimated from the observed slope of Δτ S (t) over Δt S of the target events: 15 April 2005 (light-blue line) and 2 July 2008 (red line), respectively. The arrival times of the P, S, and LP Rayleigh waves are labeled "P," "S," and "LR" ("t LR "), respectively.
constant for the entire wave packet (Figure 2d). The fractional velocity change for the LP Rayleigh wave (δV LR ) is computed as follows: where τ LR is the mean of the τ LR whose lapse time segments that satisfy the SNR and CC coefficient thresholds and t LR is the Rayleigh wave absolute arrival time ( Figure 2c). The 2σ standard deviations in δV LR are twice the computed RMS of τ i LR − τ LR ( Figure 3 and the following figures).
We emphasize that δV S and δV LR are used to highlight the time-dependent discrepancy in late 2007. If the upper crust experiences partial or uniform velocity perturbations induced by earthquake slip, then one would anticipate detectable changes in the direct S wave and in both the P and S coda waves. We observe that the τ P values are smaller than the τ S values and only detectable immediately after the Nias earthquake (see section 3). Therefore, we only focus on δV S and δV LR in our analysis. We consider the effect of fractional  Table S1). In (b), an enlarged view is displayed in the inset to better present the recovery of δV LR (open circles) during the 2005-2006 time period. The (c) δV S and (d) δV LR distributions, measured from the REs that have been relocated (Table S2) and possess high CC coefficients, and determined using Equations 5 and 6. The contrast in δV before and across the 2007 Bengkulu earthquakes, which is inferred from the probability density function for all of the data, resembles the δV distributions constructed from the individual RE sequences. changes in both the P and S wave velocities of the crust and how the velocity model perturbations influence the synthetic τ(t) results (see section 4).

Temporal Changes in δV S and δV LR : Observations of Bifurcation
We examine the τ S (t) from the N1 sequence at station PSI to determine the decade-long temporal evolution of δV S and δV LR (Figures 2 and 3). The τ S values of the first target event (15 April 2005) increase monotonically to τ S = +0.07 s (light-blue curve in Figure 2b) relative to the 9 October 2006 reference event, which is equivalent to δV S = −0.1% (dashed black line in Figure 2b). The slope of τ S decreases with time. However, the τ S values for the 2 July 2008 target event decrease to about −0.02 s (red curve in Figure 2b), which is  Similar discrepancies between τ S and τ LR are consistently detected among all of the Sumatra and Nias RE sequences at station PSI after late 2007 ( Figure 4). Furthermore, some of the P codas lag times (τ P ) and τ S measurements reveal detectable seismic-velocity changes in the P wave coda and the direct S wave onset during the first week following the 2005 Nias earthquake. Note that the τ(t) values for the S1 and S4 Sumatra RE sequences possess smaller values and fluctuate around 0 for the HF coda waves (Figures 4e and 4g). The monotonic increase in τ P only holds for the target events that occurred within the first week following the 2005 earthquake, with calculated τ P values of up to +0.03 s (thick light-blue curves in Figures 4c and 5b-5g). Because the τ P values are a factor of 2 smaller than the τ S values (e.g., Figure 3a of W. Yu et al., 2013a) to begin with, given the same decay rate, the τ P values largely decreased to zero after the first week following the 2005 mainshock (Figures 4a and 5a). For the target events that exhibit nonzero τ P , the τ S values at the direct S wave onset are of about +0.03 s (thick light-blue curves in Figures 4c and 5b-5g).
We use a kernel density estimation (Parzen, 1962) with adequate bandwidths to assess the overall distribution of the δV S and δV LR measurements whose time separations are longer than half a year and establish the consistency and robustness of the observed δV S and δV LR trends (Figures 3c and 3d). The kernel density estimation is a representation of the probability density function and can be formulated as follows: where x i is the ith δV data of the total sample size n, h is the bandwidth value that controls the smoothness of the resulting probability density curve, and K (x i , h) is the kernel expressed in Gaussian function: The peaks correspond to the maximum likelihood values for each set of δV data. Note that the break from the δV LR recovery after late 2007 is likely caused by both the 2007 Bengkulu and 2008 Simeulue earthquakes (to be discussed in section 5). We adopt the timing of the 2007 earthquakes to sort the δV values into two groups, before (blue curve, Group 1) and across the 2007 earthquakes (red curve, Group 2); the 2007 earthquakes were simply chosen for this division since they occurred first. We observe a noticeable contrast between the Group 1 and Group 2 δV values. The peak δV S value changes from −0.04% in Group 1 to +0.007% in Group 2, indicating continuous velocity recovery (Figure 3c), whereas the peak δV LR value changes from −0.16% in Group 1 to −0.25% in Group 2, indicating an additional velocity reduction and minimum recovery after 2007 (Figure 3d). The post-2007 δV LR values (Group 2) are typically a factor of 2 larger than those measured earlier in the time series (Group 1). The steady δV S recovery at Station PSI after the 2004/2005 earthquakes is consistent among all of the Sumatra and Nias RE sequences (left panels of Figures 6 and 7). The δV LR values consistently exhibit substantial offsets of −0.3% to −0.6% for all of the RE sequences across the late 2007 timeframe (right panels of Figures 6

Finite Difference Synthetics and τ(t) Investigation
The two-dimensional elastic finite difference (FD) method ) is used to compute the synthetic seismograms. Our goal is to explore the seismic velocity model parameters that are sensitive to the first-order pattern of the observed τ S and τ LR trends. Our one-dimensional (1-D) reference velocity model, termed "mST2," is a combination of the ST2 model (upper 50 km) (Lange et al., 2010) and isotropic Preliminary reference Earth model (PREM; 50-to 370-km depth) (Dziewonski & Anderson, 1981). The ST2 velocity model is derived from regional P and S wave traveltimes whose phases sampled the Sumatra Subduction Zone. Here we perturb the P wave velocity (α), S wave velocity (β), and layer thickness of the mST2 model for our sensitivity analysis (Figures 8a and 8b). Random media are adopted for all of the models, which are defined by Gaussian correlation functions with RMS perturbations and correlation lengths in the x and z directions (clx and clz, respectively) (Frankel & Clayton, 1986). A two-dimensional (2-D) random-media reference velocity

Journal of Geophysical Research: Solid Earth
model is formed from the 1-D mST2 velocity model using the following parameters: RMS = 10%, clx = 0.5 km and clz = 0.5 km (isotropic, with an aspect ratio clx/clz = 1; see Table 1). We consider two end-member models, "deep layer" (DL) and "surface" (SF). The settings for the DL models include superposed background velocity dα and dβ perturbations and variations in the clx/clz aspect ratio and RMS in the upper 1-16 km of the crust, whereas the settings for the SF models include superposed dβ perturbations, clx/clz, and RMS variations that are confined to the uppermost 0.5-1.0 km. The perturbed parameters for each model are presented in Table 1.
We consider the geometry of the N1 sequence recorded at station PSI for most of the computations, with a 22-km focal depth for the N1 RE sequence and 267-km epicentral distance from Station PSI. The dimensions of the FD regime are 500 km × 220 km (x direction × z direction), and the length of the time window is 170 s, with a 0.0625-km grid spacing and 3.125-ms time step used during the computation. The source is positioned at x = 100 km and z = 22 km. The random-media regime extends from x = 101.5 km to x = 500 km and from z = 0 km to z = 1 km for SF models and from z = 1 km to z = 16 km for DL models (Figure 8c). The FD synthetics are accurate up to 6.0 Hz. We also present two examples where we consider the geometry of the S1 sequence at Station PSI, with a 48-km focal depth for the S1 RE sequence and 534-km epicentral distance from station PSI. The dimensions of the FD regime are 900 km × 370 km (x direction × z direction) for the S1-PSI geometry, and the length of the time window is 280 s, with a 0.125-km grid spacing and 6.25ms time step. The source is positioned at x = 100 km and z = 48 km. The precision of the synthetic seismograms for S1 is accurate up to 3.0 Hz. clz is fixed at 0.5 km, meaning each scatter includes four and eight grid cells in the z direction for the S1-PSI and N1-PSI regimes, respectively. We use a dip-slip focal mechanism and 1-s triangle source time function. The HF synthetic seismograms are band-pass filtered between 0.4 and 1.6 Hz. The synthetic τ(t) values that are computed from the synthetic seismograms are based on the perturbed models relative to the reference model using a sliding-window CC approach, which is identical to the approach used to obtain the observed τ(t) values. Note that synthetic scattered wavefields are computed with the 2-D elastic FD method, and intrinsic attenuation (Q) is not included in the computation. The 2-D nature and lack of attenuation are probably the factors that contribute to the dissimilarity between the synthetic  Table 1. and observed coda waves and oscillating synthetic τ(t) in the later part of the codas. We focus on characterizing the behaviors of the HF synthetic τ S (t) bounded by the vertical dotted lines in Figures 9-12.
The heterogeneous regions that can produce the observed time-dependent time delay include a localized velocity reduction near either the source or station and a velocity reduction in the bulk of the upper crust. We first examine the effect of a localized velocity reduction near either the source or station. The N1 DL-LSRC1 model produces τ P values of up to +0.05 s and null τ S values for a localized velocity reduction near the source (orange curve in Figure 9a), whereas the N1 SF-LSRC1 model produces null τ P and τ S values (light-blue curve in Figure 9a). The N1 SF-LSTN1 model produces null τ P and τ S values for a localized velocity reduction near the station (thick green curve in Figure 9a), and the S1 SF-LSTN1 model produces oscillating τ P and τ S values between −0.05 s and +0.05 s for the S1-PSI path (thick green curve in Figure 10a). Amplifying the localized S wave velocity reduction near the station from dβ = −4% to −16% produces a linear increase in τ S (t) values, with increases of +0.075 s (N1 SF-LSTN2; thick yellow curve in Figure 9a) and +0.15 s (S1 SF-LSTN2; thick yellow curve in Figure 10a), respectively, for the N1 and S1 sequences recorded by station PSI. Two lines of argument suggest that the observed τ(t) values of the HF S coda waves are not caused by a localized near-surface velocity reduction at station PSI.
(1) The perturbed models with localized velocity perturbations near the station produce a linear increase in τ S values, even at different distance ranges, whereas the observed τ S values for the S1 sequence at station PSI does not exhibit an increase in τ S (t) (light-blue curve in Figure 4e).
(2) A monotonic increase in τ P values is detectable for the target Table 1 List

10.1029/2020JB019794
Journal of Geophysical Research: Solid Earth events that occurred within the first week following the Nias mainshock (Figures 4c and 5b-5g), but the smaller τ P values largely diminished to 0 after the first week following the Nias earthquake. The observed τ P and τ S values suggest that the velocity reduction occurs in the bulk part of the upper crust. Moreover, the localized velocity reduction near either the source or station produces τ LR values of < +0.03 s (Figures 9c and 10c).
Here we discuss how the velocity reduction parameters, RMS perturbations, and clx/clz aspect ratio for the SF and DL models influence the synthetic τ(t) values for the HF S coda and Rayleigh waves. The following FD synthetics computations employ the N1-PSI source-receiver geometry. Figure 11 displays the synthetic test results where an individual parameter is moderately perturbed, such that the sensitivity of each perturbed parameter can be assessed. Models SF1, SF2, and SF3 have a 3% RMS amplification, increase in the clx/clz ratio from 1 to 10, and additional dβ = −4% perturbation relative to the reference model, respectively. The synthetic τ(t) values for the HF S coda waves that propagated through these models all exhibit subtle fluctuations around 0, with episodic oscillations at 10, 20-30, and 47-52 s (Figure 11a). Models DL1, DL2, and DL3 have a 3% RMS amplification, increase in the clx/clz aspect ratio from one to ten, and additional dα = −0.06% and dβ = −0.12% perturbations relative to the reference model, respectively. The 3% RMS amplification also produces subtle fluctuations in the τ P and τ S values (DL1, light-blue curve in Figure 11. HF FD synthetic seismograms and corresponding τ(t) using several velocity models: (a, b) "surface" (uppermost 1.0 km, denoted as "SF") and (c, d) "deep layer" (1-to 16-km depth, denoted as "DL") with various RMS, correlation length, dα, and dβ perturbations in the random media. Synthetic LP τ LR based on various (e) SF and (f) DL model perturbations. The LP synthetics in (f) that are computed based on the reference model and model DL4 are indicated by black and green traces, respectively. The perturbed parameters of each presented model are listed in Table 1 Table 1. This figure is similar to Figure 11. 10.1029/2020JB019794 Figure 11c). However, increasing the clx/clz aspect ratio from one to ten produces a negative τ S value of −0.07 s at~50 s (DL2, orange curve in Figure 11c), and the dα = −0.06% and dβ = −0.12% perturbations produce an increase in the τ S value of up to +0.13 s (DL3, pink curve in Figure 11c). Model DL4, which includes all of the perturbations in models DL1, DL2, and DL3, produces a monotonic increase in τ P , as well as smaller τ S values at~50 s (green curve in Figure 11c). The smaller τ S values mimic the negative τ S values in the DL2 model at 47-52 s due to the increase in the clx/clz aspect ratio.
We now consider the effects of RMS amplification and perturbing several parameters within a given velocity model. A monotonic increase in both τ P (t) and τ S (t) is still better approximated by the DL models for the HF P and S wave codas (DL5 and DL6 in Figure 12c). The SF models with a strong S wave velocity reduction (dβ = −16%) produce oscillating τ P (t) and τ S (t) values with peaks at 10, 20-30, and 47-52 s (SF5 and SF6 in Figure 12a). The DL models with RMS = 25%, dα = −0.06%, dβ = −0.12% perturbations, and clx/clz = 1 produces linear increases of up to +0.072 and +0.147 s in τ P and τ S , respectively (DL5, light-blue curve in Figure 12c; the 2D dβ profile is displayed in Figure 8d). The τ P values increase by a factor of 1.5 when the clx/clz aspect ratio is increased to 10, with a value of up to 0.11 s; the τ S values are larger up to +0.138 s over the 40-to 45-s window, but then decrease over the 47-to 52-s window (DL6, orange curve in Figure 12c; the 2-D dβ profile is displayed in Figure 8e). This τ S pattern of an increase followed by a drop over the 47-to 52-s window is similar to that observed for models DL3 and DL4 (Figure 11c). The heterogeneous region in model DL7 only spans 50% of that in model DL5; its predicted τ P values exhibit an increase over the 10-to 20-s window to +0.025 s, followed by a drop over the 20-to 30-s window to −0.03 s, whereas the τ S values exhibit a monotonic increase up to +0.1 s (DL7, pink curve in Figure 12c). The combined perturbations of the SF and DL models produce τ P and τ S values that are similar to those based calculated from the DL models, with larger τ P and smaller τ S values often obtained ( Figure 12e; the 2-D dβ profile of model SF4 + DL7 is displayed in Figure 8f). The slope of synthetic τ S (t) based on the above DL and combined perturbations of the SF and DL models estimate δV S in the range of −0.61% (DL6) and −0.98% (DL7), slightly lower than the observed δV S immediately after the 2005 event (Figures 3 and 6). Models DL5 and SF4 + DL5 produce δV S = −0.72% and −0.86%, respectively, indicated as the slope of dotted line in Figures 12c and 12e.
Our modeling results indicate that a strong S wave velocity reduction in the near surface substantially contributes to the large τ LR values but has a minimal effect on the observed τ P and τ S patterns. A subtle S wave velocity reduction in the upper crust is necessary to account for the observed τ P and τ S patterns and can also partially contribute to the τ LR pattern. Depending on the availability of early τ S and τ LR measurements after the 2004/2005 mainshocks, the amount of the S wave velocity reduction in the near surface and upper crust can be scaled accordingly. Our preferred explanation is as follows. The observed τ S and τ LR patterns after the 2004/2005 earthquakes may reflect the combination of a strong S wave velocity reduction in the near surface and a subtle S wave velocity reduction in the upper crust. However, the τ LR amplification and the null τ S after late 2007 predominantly reflects another substantial S wave velocity reduction in the near surface (red curves in Figure 4). The scattering media changing from isotropic (clx/clz = 1) to vertically transverse isotropic (clx/clz = 10) and amplifying the RMS can influence S coda and Rayleigh wave propagation, thereby altering the τ S and τ LR values. We discuss the likely scenarios that can potentially explain these observations and modeling results in the following section.

δV S Recovery After the 2004 and 2005 Earthquakes
Obermann et al. (2013) positioned the source and receiver at the surface in their coda wave computations and suggested that the later part of the coda waves were predominantly sensitive to the bulk velocity heterogeneities of the media, with the sensitivity of these waves to the source depth being approximately bounded by the length scale of the mean free path. Here we position the source and receiver at 22-km depth and the surface, respectively, in our FD simulations. Our FD synthetic tests reveal that the DL models with subtle S wave velocity perturbations that are confined to the upper 1-16 km produce characteristics that are consistent with a monotonic increase in τ S (Figure 12c), whereas the SF models with near-surface velocity perturbations (confined to the uppermost 0.5-1.0 km) cannot reproduce this key feature (Figure 12a). Although the source-receiver geometries of the two modeling approaches are different, our inferences on the sensitivity of the later part of the coda waves to the bulk heterogeneities of the media are consistent with those in Obermann et al. (2013).
The δV S observations show continuous recovery from 2005 to 2015. It is unlikely that the δV S induced by coseismic slip would be sustained for several years, even though the deformation induced by large earthquake coseismic slip can produce structural changes in the bulk crust and time-varying δV S behavior. Additional independent arguments suggest that the continuous δV S recovery from 2005 Table S2). The postseismic displacement time series recorded by GPS Stations UMLH and LEWK, and the GPS normalized displacement time series relative to the concurrent δV S values are provided in Figures S4 and S5, respectively, for comparison. Note that the details of the GPS data processing can be found in Feng et al. (2015). The logarithmic recovery of δV S corresponds to the postseismic deformation time series reasonably well. Feng et al. (2015) suggested that viscoelastic relaxation was essential for explaining the long-term geodetic measurements. Qiu et al. (2018) recently estimated a background upper mantle viscosity of 10 17 -10 19 Pa s beneath the Sumatra region, which would indicate relaxation timescales (days to years) that are similar to our observed timescales of δV S recovery.

Multiple δV LR Reductions Induced by the 2004, 2005, 2007, and 2008 Earthquakes
The observed δV LR values at Station PSI from late 2007 to early 2008 exhibit a remarkable reduction of at least −0.3%, which is a factor of 2 larger than the reductions caused by the 2004 and 2005 earthquakes (Figures 3, 6, and 7). These new data reveal that the δV LR (Hobiger et al., 2016;Rubinstein & Beroza, 2004b). We examine the PGV of the large earthquakes recorded by Station PSI in the 0.5-to 2.0-Hz frequency band (Table S3). The PGV amplitude of the 2005 M w 8.6 earthquake is larger than that of the 2004 M w 9.2 event, and the PGV amplitudes of the 2007 M w 8.4 and M w 7.9, and 2008 M w 7.3 earthquakes are smaller than those of the 2005 M w 8.6 event, even though a stronger δV LR reduction that is associated with the 2007 and 2008 earthquakes is detected. We speculate that the strong GM induced by the 2004 and 2005 events made the subsurface cracks and fractures more susceptible to the GM induced by the 2007 and 2008 earthquakes, resulting in another substantial velocity reduction.

Temporal Changes in δV During 2005-2015
REs discretely sample δV with a limited temporal resolution. However, concurrent δV S and δV LR observations, and the differences in their temporal patterns offer new insights into the diverse processes that modulate media properties. Figure 13 summarizes the decadal δV S and δV LR observations and illustrates our inferences on the evolution of the crustal seismic-velocity structure along the Banda Aceh-Nias-Simeulue segment of the Sumatra Subduction Zone. A subtle −0.1% δV S reduction due to the 2004 and/or 2005 earthquakes is followed by steady recovery, with no noticeable disruption from 2005 to 2015 (δV S in Figure 13). However, a −0.1% δV LR reduction is induced by the 2004 earthquake, and a slightly larger −0.1 to −0.3% reduction is induced by the 2005 event. We find that δV LR experiences recovery until late 2007, with the GM from the 2007 and 2008 earthquakes then inducing a substantial δV LR reduction (−0.3% to −0.6%) that possibly occurred in multiple steps. We speculate that any δV LR caused by the 2010 Banyak earthquake may have gone undetected due to the limited temporal resolution of the data, resulting in the observed gentle velocity recovery trend.
In summary, the differences in the temporal evolution of δV S and δV LR inferred from REs are a powerful tool for discriminating between the mechanisms that drive seismic-velocity changes after great earthquakes. Steady δV S recovery is consistent with the logarithmic time dependence inferred from postseismic GPS data ( Figures S4 and S5), which predominantly reflects afterslip that occurred in the updip and downdip sections of the rupture area (Feng et al., 2015) ( Figure 1) and/or viscoelastic relaxation after the 2004 and 2005 earthquakes. Furthermore, two-dimensional FD synthetics computations demonstrate that the monotonic increase in the observed τ S values can be explained by a subtle velocity reduction of dβ = −0.12% at 1-to 16-km depth (DL5 in Figure 12c). However, δV LR displays additional drops after the large earthquakes during the 2007-2008 time period, breaking the pattern of monotonic recovery. The model with a subtle S wave velocity reduction at 1-to 16-km depth can partly account for the observed τ LR pattern after the 2004/2005 events, but it is insufficient to produce values that are comparable to the observed τ LR pattern after late 2007. Strong S wave velocity reductions of dβ = −4% and −16% that are confined to the uppermost 0.5-1.0 km can generally account for the observed τ LR values after 2005 and late 2007, respectively (Figure 12g). These sudden changes are primarily associated with a cycle of damage, healing, and redamage in the subsurface layer (Sawazaki et al., 2015(Sawazaki et al., , 2018Vidale & Li, 2003).

Data Availability Statement
The seismic data analyzed in this paper were assembled from the PS network, which were downloaded from the JAMSTEC data center, Japan (http://p21.jamstec.go.jp/top/; last accessed on 23 August 2019).