Seismic and Aseismic Fault Slip During the Initiation Phase of the 2017 M W = 6.9 Valparaíso Earthquake.

Transient deformation associated with foreshocks activity has been observed before large earthquakes, suggesting the occurrence of a detectable preseismic slow slip during the initiation phase. A critical issue consists in discriminating the relative contributions from seismic and aseismic fault slip during the preparation phase of large earthquakes. We focus on the April-May 2017 Valparaíso earthquake sequence, which involved a M W  = 6.9 earthquake preceded by intense foreshock activity. To assess the relative contribution of seismic and aseismic slip, we compare surface displacement predicted from foreshocks source models with transient motion measured prior to the mainshock. The comparison between observed and predicted displacements shows that only half of the total displacement can be explained by the contribution of foreshocks. This result suggests the presence of aseismic pre-slip during an initiation phase preceding the mainshock.

Introduction

Experimental and theoretical studies suggest that earthquakes begin with aseismic slow slip accelerating into a dynamic, catastrophic rupture (Das & Scholz, 1981; Kaneko et al., 2016; Latour et al., 2013; Ohnaka, 2000). Laboratory‐derived rate‐ and‐state models depict different evolution of pre‐slip within nucleation zones of various sizes (Ampuero & Rubin, 2008; Kaneko & Ampuero, 2011). With technological advances such as high‐speed photoelastic techniques, the progressive acceleration from slow stable slip to fast dynamic slip can be accurately monitored in laboratory conditions (e.g., Latour et al., 2013). Despite these advances, the detectability of such nucleation phases on natural faults is still an open question. In addition to the nucleation itself, observations of the precursory phase leading to an earthquake indicate that earthquakes are often preceded by foreshocks that could potentially be triggered by aseismic pre‐slip (Bouchon et al., 2011, 2013; Kato et al., 2012). Nonetheless, the role of foreshocks during this precursory phase remains unclear. At present, two end‐member conceptual models compete in explaining the occurrence of foreshocks. In the first model, foreshock stress changes contribute to a slow cascade of random failures, leading eventually to the mainshock (Ellsworth & Bulut, 2018; Helmstetter & Sornette, 2003; Marsan & Enescu, 2012). The second model proposes that foreshocks are triggered by aseismic slip corresponding to the nucleation process of the mainshock (Bouchon et al., 2011; Dodge et al., 1996).

The continued development of geophysical networks in active tectonic regions provides new opportunities to better capture the genesis of earthquakes. Geodetic observations provide strong evidences of preseismic transient deformations at various time‐scales (Ito et al., 2013; Mavrommatis et al., 2014; Ozawa et al., 2012; Socquet et al., 2017; Yokota & Koketsu, 2015). However, the interpretation of such observations is often difficult. This is particularly evident for the 2014 M  = 8.4 Iquique (Chile) earthquake, which was preceded by an active foreshock sequence that started 8 months before the mainshock (Kato & Nakagawa, 2014). This foreshock sequence was accompanied by clear Global Positioning System (GPS) transient displacements, corresponding at least to some extent to aseismic fault slip preceding the mainshock (S. Ruiz et al., 2014; Socquet et al., 2017). The aseismic behavior of the observed preseismic transient is however debated as it might largely correspond to the cumulative co‐seismic displacement of the foreshocks and associated afterslip (Bedford et al., 2015; Schurr et al., 2014). A reliable estimate of the relative contribution of seismic and aseismic deformations during nucleation is essential to better capture fault processes at the onset of earthquakes (Herman et al., 2016).

On April 24, 2017, a M  = 6.9 earthquake occurred offshore Valparaíso in the central segment of the Chilean megathrust (33.089°S, 72.116°W, 21:38:28 UTC; Centro Sismólogico National, CSN). This event is relatively moderate given that this region of the Chilean subduction experienced earthquakes of magnitudes M  > 8 (Comte et al., 1986; Dura et al., 2015). This earthquake, however, caught the attention of seismologists because it was preceded by a vigorous foreshock activity in the ∼2 days preceding the mainshock. This precursory activity has also been captured by GPS stations indicating a preseismic trenchward motion over a similar time‐scale (S. Ruiz et al., 2017; J. A. Ruiz et al., 2018). A preliminary analysis of seismological and geodetic observations suggests that 80% of preseismic GPS displacement is due to aseismic fault slip preceding the mainshock (S. Ruiz et al., 2017). This first order estimate is obtained by comparing inverted pre‐slip with the seismic moment of foreshocks assuming they are all located on the subduction interface. This assumption is questionable as seismicity catalogs depict a significant dispersion of earthquake locations around the plate interface (S. Ruiz et al., 2017; J. A. Ruiz et al., 2018), most events being located at depths larger than the slab 1.0 model (Hayes et al., 2012). Such dispersion, probably related to depth uncertainty, implies a significant nonrandom bias in seismic moment for dip‐slip earthquakes. For example, if an earthquake at 20 km depth is mislocated at 25 km, the moment is underestimated by nearly 20% using long‐period teleseismic records (Tsai et al., 2011). Such mis‐estimation of seismic moment may lead to nonnegligible errors in the contribution of foreshocks to observed preseismic deformations.

The primary goal of this study is to assess the relative contribution of seismic and aseismic slip during the few days preceding the 2017 Valparaíso earthquake. Estimating the seismic contribution to observed geodetic displacement is difficult as we deal with moderate‐sized foreshocks (M  < 6) for which a co‐seismic offset is not clearly visible on GPS time‐series. The seismic contribution to the observed displacement can be estimated by modeling the source of foreshocks from seismic data. However, this process should be done carefully as source models and the corresponding predictions can be affected by significant uncertainties. In this study, we obtain a moment‐tensor catalog and predict the corresponding co‐seismic offsets at GPS stations accounting for observational and modeling uncertainties. In particular, we account for prediction uncertainties associated with inaccuracies in the Earth model. We find that about half of the observed GPS preseismic displacement is aseismic and is caused by pre‐slip in the vicinity of the impending mainshock hypocenter. Such preseismic deformation is unlikely to be explained by afterslip induced by preceding foreshocks. This suggests that aseismic pre‐slip played an important role in the 2017 Valparaíso sequence.

Preseismic Transient Displacements Captured by GPS

We process GPS data of 68 stations in South America from several networks (CSN, LIA Montessus de Ballore, Ministerio de Bienes Nacionales, RAMSAC, RBMC‐IP, IGS, IGM Bolivia, see supplementary information for references). Processing is done using a differential approach (Herring et al., 2018) including tropospheric delays and horizontal gradients. The results are computed in the ITRF 2014 reference frame (Altamimi et al., 2016) and converted in a fixed South‐America frame (Nocquet et al., 2014). We use daily solutions except for the last position before the mainshock, which is obtained from data up to 1 h before the event. We remove a trend corresponding to interseismic motion from the time‐series by fitting a linear regression in a 4 months time‐window before the mainshock. Finally, we subtract the first sample of the time‐series (i.e., which we consider as displacement zero) and obtain the corresponding offsets.

Figures 1b and S3 show the resulting horizontal displacements for stations in the vicinity of the study area. There is a clear westward motion starting about 3 days before the mainshock and reaching ∼8 mm close to the coast. Figure 1b compares GPS time‐series with the cumulative number of earthquakes in the micro‐seismicity catalog obtained by S. Ruiz et al. (2017). Interestingly, the preseismic GPS transient starts before a noticeable increase in seismicity. In Figure 1b, we can see that the slope of cumulative seismicity rate does not change significantly at the beginning of the transient. The increase in seismicity rate is delayed by about 24 h and only starts with a M  = 6.0 foreshock on April 23 (purple star in Figure 1b). This suggests that aseismic pre‐slip initiated on the fault before the increase in foreshock activity.

Figure 1

The 2017 Valparaíso earthquake sequence. (a) Earthquake locations including foreshocks (blue circles), mainshock (green star), and aftershocks (white circles). The red colormap indicates the pre‐slip distribution resulting from the inversion of Global Positioning System (GPS) data (see section 5). The black arrows show the cumulative observed GPS surface displacements (up to 1 h before the mainshock). Orange dots indicate the seismicity distribution from January 01, 2017 until October 05, 2017 according to the microseismicity catalog obtained by S. Ruiz et al. (2017). (b) GPS Time‐series in the vicinity of Valparaíso. The vertical red dashed line indicates approximate onset of the transient displacement visible on the time‐series. The cumulative number of earthquakes from S. Ruiz et al. (2017) is shown at the bottom of the figure. The purple star represents the largest M  = 6.0 foreshock.

Centroid Moment Tensor Catalog

To constrain the contribution of foreshocks to the observed GPS displacement, we estimate Centroid Moment Tensor (CMT) parameters for moderate to large earthquakes during the Valparaíso earthquake sequence (from April 05, 2017 up to May 30, 2017). We use records from broadband seismic stations located within 12° from the mainshock hypocenter. These stations are mostly included in the C and C1 regional networks maintained by the Centro Sismológico Nacional (CSN) of the Universidad de Chile (Universidad de Chile, 2013). We also use stations operated by GEOSCOPE, and IRIS/USGS network (Albuquerque Seismological Laboratory (ASL)/USGS, 1993, 1988; Institut de Physique du Globe de Paris and Ecole et Observatoire des Sciences de la Terre de Strasbourg (EOST), 1982).

We use a modified version of the W‐phase algorithm adapted to regional distances and the magnitude range of the Valparaíso sequence (Kanamori & Rivera, 2008; Zhao et al., 2017). Estimated parameters are the deviatoric moment tensor, the centroid location, the centroid time, and the half‐duration of an isosceles triangular moment rate function. The inversion is performed by fitting full waveforms in a 180 s time‐window starting at the P‐wave. We filter data between 12 s and 100 s using different pass‐bands for different magnitude events (see Table S1). We compute Green’s functions for the source inversion in a one dimensional layered structure extracted from the three dimensional Earth model of S. Ruiz et al. (2017) in the area of Valparaíso (Figure S4).

The resulting CMT catalog is shown in Figure 2 and in Table S2. Most earthquakes (more than 90% of the total catalog) have thrust mechanisms. Interestingly, foreshocks are mostly concentrated close to the mainshock hypocenter (see Figures 1 and 2a). On the other hand, aftershocks show a different behavior, surrounding the region where foreshocks have previously occurred.

Figure 2

Centroid Moment Tensor (CMT) solutions of the 2017 Valparaíso earthquake sequence and cumulative moment (a) CMT solutions of the 2017 Valparaíso earthquake sequence. Focal mechanisms are contoured in blue and black for foreshocks and aftershocks respectively. The size of beach balls scales with the moment magnitude. Color of the compressive quadrants represents the event depth. (b) Cumulative scalar seismic moment of the 2017 Valparaíso sequence. The mainshock scalar moment is not included in this figure. The red dashed line outlines the approximate onset of transient displacements visible on GPS time‐series. The green line indicates the origin time.

The cumulative scalar seismic moment released by foreshocks before the mainshock is largely dominated by two events with M  ≥ 5.5 (cf., Figure 2b). These foreshocks of magnitude M  = 6.0 and M  = 5.5 occurred respectively 43 h and 26 h before the mainshock. As our CMT catalog only consists of M  ≥ 3.8 earthquakes, the contribution of microseismicity is not included in our estimates of cumulative seismic moment before the mainshock. Even though the individual contribution of these small earthquakes to the observed displacement is negligible, their large number may contribute to surface displacement. To assess the contribution of small earthquakes, we consider the frequency‐magnitude distribution of our CMT catalog assuming a completeness magnitude of M  = 3.9 (Figure S5). We compare our catalog with previous moment tensor catalogs of the same sequence (S. Ruiz et al., 2017; J. A. Ruiz et al., 2018), which are qualitatively consistent with our estimates (Figure S5). We then compute the Gutenberg‐Richter (GR) law using the methodology proposed by Aki (1965) for the whole sequence, and the foreshocks sequence. Even though the GR laws show some discrepancies, they are in good agreement considering the uncertainties on our estimates (Figure S5). The foreshocks GR law is then extrapolated to lower magnitudes, and the cumulative moment of magnitudes below the magnitude of completude is included to correct for the influence of small, hence not detected earthquakes. Our CMT catalog suggests a cumulative moment M 0 = 1.474 × 1018 N⋅m. The cumulative seismic moment of foreshocks with magnitudes below completeness is M 0 = 4.966 × 1015 N⋅m (i.e., M  = 4.4). The contribution of microearthquakes is therefore negligible compared to seismic events.

To evaluate the contribution of foreshocks to observed surface displacements, we calculate synthetic static displacements using our CMT catalog and the same one dimensional velocity model employed to obtain our CMT solutions. Synthetics are computed using the Classic Slip Inversion (CSI) package (http://www.geologie.ens.fr/jolivet/csi) incorporating the approach of Zhu and Rivera (2002) to compute static displacement in a layered model. Results on Figure S6 indicate that the largest foreshock (M  = 6.0) largely dominates the co‐seismic contribution to the observed GPS transient while M  < 6.0 events in our catalog generate relatively small surface displacement. Assuming that microearthquakes are located in the vicinity of M  ≥ 3.8 foreshocks, they should also have a negligible contribution to the observed surface displacement (given their small cumulative scalar moment). As the M  = 6.0 foreshock plays a important role in the sequence, we assess uncertainties associated with the corresponding CMT parameters.

Uncertainty on Predicted Co‐seismic Displacements

Synthetic co‐seismic surface displacements are sensitive to uncertain earthquake source parameters. For large magnitude foreshocks, uncertainties on centroid location and moment tensor affect our estimates of the co‐seismic contribution to the transient displacement observed before the mainshock. Source parameters uncertainties can either result from observational errors, or from errors in the forward model (prediction/theoretical errors). For example, there might be innaccuracies in the velocity model, which is known to induce nonnegligible errors in CMT solutions (Duputel et al., 2012, 2014; Morales‐Yañez et al., 2020). The point source assumption is another source of uncertainty in the forward model. As for the observations, temporally and spatially variable noise level at seismic stations is a major source of uncertainty.

In order to assess uncertainties associated with the CMT solution of the largest M  = 6.0 foreshock, we perform a new CMT inversion within a Bayesian framework, following Duputel et al. (2012, 2014). Each source of uncertainty considered here is integrated in the problem as a covariance matrix. The covariance matrix C , associated with observational errors, is derived after a first CMT inversion. From this inversion, an average correlation function is derived from residuals between synthetic and observed waveforms at each station. This allows us to estimate the correlation between neighbor data samples, and include it into C . The standard deviation for each channel is fixed to 4 times the corresponding average absolute residuals. This empirical procedure provides a conservative estimate of observational uncertainty associated with each waveform.

Forward modeling uncertainties are represented by the matrix C , which assesses the influence of inaccuracies in the Earth model. We use the same velocity model as in section 3 assuming log‐normal uncertainties on elastic parameters as shown in Figure S4. Uncertainty in each layer is estimated by assessing the spatial variability of the three dimensional Earth model of S. Ruiz et al. (2017) in the epicentral region and by comparison with other regional models (e.g., J. A. Ruiz et al., 2018). To evaluate the corresponding variability in the predictions, we employ the first‐order perturbation approach described in Duputel et al. (2014), assuming that prediction error is linearly related with uncertainty on the elastic parameters. A test is described in supplementary information S2 and Figures S7 and S8 to assess the validity of this approach.

The posterior ensemble of plausible source locations and moment tensors is appraised using a strategy similar to Sambridge (1999). At a fixed point‐source location in time and space, the posterior distribution of moment tensor parameters is Gaussian and can be written as (Tarantola & Valette, 1982): where m are the moment tensor parameters, d is the data vector containing the concatenated observed waveforms and x is the point source location. The right‐hand member of this equation is a Gaussian distribution of mean and covariance . The posterior mean is the maximum a posteriori moment tensor given by: where G is the Green’s function matrix while C  = C  + C is the covariance matrix reflecting observational (C ) and prediction uncertainties (C ). The posterior covariance matrix is given by:

To get the joint posterior distribution on moment tensor m and source location x, we first calculate and on a three dimensional grid of possible point‐source locations around the hypocenter. Starting from the initial location x determined in Section 3 (corresponding a moment tensor m ), we then employ an hybrid metropolis algorithm by repeating the following iterations until a sufficiently large number of model samples is generated:

Randomly generate a candidate point‐source location x* = x  + δ x where δ x is a small perturbation randomly generated from a Gaussian distribution with a standard deviation of 0.1° in latitude/longitude and σ = 0.1 km in depth

Extract and from the grid point closest to x* and generate a random model m* from p (m|d , x*) in Equation 1

Accept or reject m* and x* using a standard Metropolis approach

•Draw a random number α ∼ U (0, 1)

•Accept m* and x* if

•Otherwise duplicate m and x

Figure 3 shows 4,500 model samples generated using the approach described above. The posterior distribution shows a location uncertainty of about 10 km. We observe a good fit between observed and synthetic seismograms (Figure S9). However, we also notice a trade‐off between longitude and depth, which probably results from the distribution of stations used for inversion (Figure S10). To evaluate the uncertainty on the predicted co‐seismic displacement, we simulate static displacement for each model samples shown in Figure 3. The resulting stochastic co‐seismic displacements are shown in Gray in Figure 4a for GPS stations that are closest to the mainshock epicenter. This shows prediction uncertainties ranging from 0.25 to 0.4 mm on the east component of displacement. Despite these uncertainties, the predicted cumulative co‐seismic offsets are still significantly smaller than the observed preseismic displacements (∼ 6–8 mm of the east component for the closest stations).

Figure 3

Bayesian point‐source model for the M  = 6.0 foreshock on April 23, 2017. Blue circles and lines in the figure represent model samples randomly drawn from the posterior distribution. (a) Samples from the posterior PDF depecting uncertainties in the point source location. The red and orange stars are the initial solution (i.e., starting model) and the posterior mean model respectively. (b) Focal mechanism uncertainty. (c) Marginal posterior PDF of the scalar seismic moment. The red and orange lines are the initial and the posterior mean model.

Figure 4

Slip during the Valparaíso foreshock sequence. (a) Time series of global positioning system data (blue) and stochastic foreshock‐induced co‐seismic displacement (gray). Red dots represent the average of stochastic co‐seismic offsets. Green cross corresponds to the total foreshock displacement, including the contribution of earthquakes below the magnitude of completeness. (b) Distributions of observed preseismic displacement and predicted cumulative co‐seismic offsets caused by foreshocks. Blue histograms represent observations assuming Gaussian uncertainties from standard errors estimated at each station. Red histograms correspond to the posterior distribution of cumulative foreshock‐induced co‐seismic displacement. (c) Percentage of aseismic displacement for each station. (d) Average postseismic signal measured on stations TRPD, VALN, BN05 and QTAY (see Figure S11). (e) Slip inversion of preseismic GPS data. (f) Slip inversion of GPS data after removing foreshock‐induced displacement. Black and blue arrows are observed and predicted horizontal GPS displacements along with their 1‐σ ellipses (representing observational and prediction uncertainties, respectively). Colored circles are observed (outer circles) and predicted (inner circles) vertical displacements from GPS and tide gauges, respectively.

Partitioning Between Seismic and Aseismic Fault Slip

In Figure 4, we compare the total cumulative foreshock co‐seismic offset with the observed preseismic GPS displacement. Predicted co‐seismic displacements include the contribution of microearthquakes below the magnitude of completeness, assuming a total scalar moment derived from our GR analysis with a location and mechanism similar to the M  = 6.0 foreshock. As discussed earlier, only the largest foreshock M  = 6.0 is significantly contributing to co‐seismic displacements (see Figures 4a and S6). The contribution of earthquakes smaller than M  = 6.0 has a minimal impact on the final result.

To get a total budget of seismic and aseismic displacement before the mainshock, Figure 4b compares GPS data 1 h before the mainshock with the corresponding cumulative foreshock displacement. Observed displacement are on average between 4 and 6 mm larger than co‐seismic offsets. Such differences cannot be explained by uncertainties on the observations and the predictions. These results clearly suggest that a significant portion of the observed preseismic deformation is actually aseismic and cannot be caused by foreshocks. We estimate that about 51 ± 11% of the displacement measured at the GPS stations originates from aseismic slip on the megathrust. As shown in Figure 4c, the portion of aseismic deformation is quite consistent between stations suggesting that a common source located in the vicinity of the foreshocks could explain those results.

To further explore this hypothesis, we then conduct two inversions: a first slip inversion of the total GPS preseismic displacement and another inversion after removing the contribution of foreshocks (i.e., aseismic displacement only). To build a fault geometry, we use the CSI package to mesh the Slab 2.0 model with triangles of variable sizes as shown in Figures 4e and 4f. We invert for slip values at the triangular nodes using AlTar, a Markov chain Monte Carlo sampler based on the algorithm described by Minson et al. (2013). Continuous fault slip distribution is represented as a linear interpolation of the slip values at the triangular nodes. Green’s functions are computed in the same stratified elastic model used for our CMT catalog (Figure S4). Given the limited amount of available observations, we enforce a positive Laplacian prior distribution with a scale parameter of 1 m. Such sparsity‐inducing prior will favor ”simple” models with slip only where it is requested by the data. Results in Figures 4e and 4f shows that GPS observations can be explained by slip in the vicinity of the mainshock hypocenter. Aseismic slip distribution appears to be somewhat more spread out, which may be an effect of the larger uncertainty associated with GPS data after removing the contribution of foreshocks (as the co‐seismic prediction uncertainty propagates in the corrected GPS data).

Discussion and Conclusion

We investigate the seismic and aseismic motions during the preparation phase of the 2017 M  = 6.9 Valparaíso earthquake. We first evaluate the contribution of foreshock‐induced displacement to preseismic GPS observations. Co‐seismic offsets are largely dominated by a M  = 6.0 foreshock that occurred ∼43 h before the mainshock. As pointed out in Section 2, the transient GPS signal starts before the increase in seismicity rate. More specifically, we can see in Figure 4a that the observed displacement on April 22 mainly corresponds to aseismic slip as no significant foreshock occurs on that day. On the other hand, the position on April 23 results from a combination of seismic and aseismic fault slip. The detailed evolution of the partitioning between seismic and aseismic slip is difficult to interpret using daily GPS time‐series in which each position corresponds to an average over 24 h. This analysis is also subject to large observational and prediction uncertainties. For these reasons, we focus on the overall partitioning between seismic and aseismic slip during the preparation phase of the Valparaíso earthquake.

Our analysis shows that a significant part of preseismic GPS observations are not explained by foreshock‐induced displacement even when accounting for prediction and observation uncertainties. We estimate that ∼ 50 ± 11% of GPS displacements is likely caused by aseismic slip, a ratio that is fairly consistent for different stations in the vicinity of the Valparaíso sequence (Figure 4c). To check whether such preseismic motion could be explained by slip on the plate interface, we conduct a slip inversion after correcting GPS data from foreshock‐induced displacement (cf., Figure 4f). The distribution of aseismic preslip spreads toward the west of Valparaíso city with an extension of about 50 × 90 km and a scalar moment of M 0 = 3.08 × 1018 N.m (i.e., M  = 6.26). This aseismic motion represents about 50% of the moment calculated for the slip model derived from uncorrected GPS data (M 0 = 5.67 × 1018 N.m, Figure 4e). Given the cumulative moment of foreshocks (M 0 = 1.48 × 1018 N.m), we estimate that nearly 70% of the scalar moment released during the preparation phase of the Valparaíso mainshock is aseismic, which is roughly in agreement with estimates from S. Ruiz et al. (2017). The smaller portion of aseismic moment derived from the comparison of slip models in Figures 4e and 4f likely results from the simplistic assumption in Figure 4e that all foreshocks are located on the plate interface.

Even if our analysis demonstrates the existence of aseismic slip prior to the Valparaíso mainshock, such aseismic motion may include afterslip from preceding bursts of seismicity. This has been suggested for preseismic displacement observed before the 2014 M  = 8.1 Iquique earthquake, which could potentially be explained by afterslip induced by foreshock seismicity (Bedford et al., 2015). Testing such possibility for the 2017 Valparaíso sequence is difficult as we cannot easily isolate the afterslip signal from GPS time‐series, which likely incorporate other contributions including pre‐slip of the impeding mainshock. To assess the contribution of afterslip, we employ two approaches. In a first approach, we use the mainshock postseismic GPS signals as a proxy for the afterslip induced by foreshocks. The mainshock postseismic time‐series are normalized by the co‐seismic offset of each station to evaluate the relative proportion of postseismic displacement as a function of time. This suggests that about 10% of the co‐seismic moment after 43 h corresponds to postseismic deformations (see Figures 4d and S11). This result is consistent with values reported for earthquakes with similar or larger magnitudes (Chlieh et al., 2007; D’agostino et al., 2012; Lin et al., 2013). If we assume a similar behavior for the foreshocks, the postseismic signal caused by foreshocks is below measurement uncertainties (∼0.7 mm for an uncertainty of 1.1 mm in GPS signals) and can therefore be neglected. In a second approach, we make the more conservative assumption that afterslip caused by foreshocks is totally released before the mainshock. Following the empirical scaling relationship M 0(postseismic)/M 0(coseismic) = 0.36 ± 0.2 proposed by Alwahedi and Hawthorne (2019), the aseismic displacement not related to foreshocks is reduced to about 37% ± 13% of the total preseismic GPS observations (Figure S12). The total observed displacement is therefore unlikely to be explained by the contribution of foreshocks even when adding the associated afterslip. Such evaluation should be taken with caution due to the nonlinear nature of the relationship between slip rate and co‐seismic stress change for afterslip (e.g., Perfettini & Avouac, 2004; Perfettini et al., 2010).

Diverse numerical and experimental studies bring up the potential importance of aseismic pre‐slip in the triggering of foreshocks (e.g., Kaneko et al., 2016; McLaskey & Kilgore, 2013). If such observations apply on natural faults, foreshock locations could potentially inform us about the overall spatial extent of the nucleation zone prior to an earthquake. This idea is in fairly good agreement with our results suggesting a first‐order correlation between pre‐slip distribution and the location of foreshocks (Figures 1 and 4). Even if pre‐slip appears to be an important mechanism in the triggering of foreshocks, part of the foreshock activity likely results from cascading phenomena due to stress changes of neighboring events. In addition, we still need to understand why most earthquakes are not preceded by foreshock activity and even less with observable preseismic motion. This lack of systematic precursory activity might partly be due to an observational gap due to the incompleteness of current seismicity catalog (as suggested by Mignan, 2014) or the lack of near fault geodetic observations prior to large earthquakes. The analysis of an highly complete earthquake catalog in Southern California showed that 72% of M  ≥ 4 earthquakes in the region are preceded by an elevated seismic activity compared with the background seismicity rate (Trugman & Ross, 2019), suggesting that foreshock activity is more ubiquitous than previously thought. However, a recent reanalysis of the same catalog suggested that a much smaller portion of these foreshock sequences were really anomalous and could not be attributed to temporal fluctuations in background seismicity rate (van den Ende & Ampuero, 2020). Although anomalous foreshock sequences currently appear to be the exceptional, the improvement of near‐fault geodetic and seismological observational capabilities are essential to bridge the gap between natural fault observations and laboratory experiments, where foreshocks are commonly observed.

Supporting Information S1

Click here for additional data file.