Megathrust earthquakes induce local and large-scale deformations (e.g., Hu et al., 2004; Vigny et al., 2011; Moreno et al., 2012), which may last several decades due to viscoelastic relaxation (Khazaradze et al., 2002; Hu et al., 2004; Wang et al., 2012). The crust continues to move after the coseismic instantaneous response, exhibiting different deformation patterns that are observed by GPS measurements (Wang et al., 2012; Tomita et al., 2015; Klein et al., 2016). Inland GPS data show that the crust moves in the sense of the main shock in front of the rupture area, which in most subduction zones represents trenchward movement. At neighboring segments of the rupture area, however, the GPS horizontal displacements exhibit landward motion (Heki and Mitsui, 2013; Tomita et al., 2015; Klein et al., 2016; Loveless, 2017; Melnick et al., 2017), as envisioned by Anderson (1975). While the trenchward motion patterns are mostly attributed to afterslip on the fault interface and viscous relaxation in the lower crust and upper mantle (e.g., Hu et al., 2004; Hergert and Heidbach, 2006; Wang et al., 2012; Peña et al., 2019, 2020), the driving mechanisms controlling landward patterns following megathrust earthquakes are currently a matter of controversial debate. For instance, Hu et al. (2016) and Tomita et al. (2015) showed that landward acceleration of the oceanic plate following the 2011 Tohoku-oki earthquake in Japan may be attributed to upper mantle viscoelastic relaxation. Alternatively, other authors have proposed that the slab pull balance of forces after large earthquakes is the driving mechanism of the GPS-observed postseismic landward patterns (Heki and Mitsui, 2013; Yuzariyadi and Heki, 2020).

The deformation patterns induced by the coseismic and postseismic deformation are of vital importance, as they can potentially modify the state of stress at adjacent segments, which may be mature segments critically stressed, and ultimately triggering large-magnitude earthquakes. Examples of earthquake triggering because of coseismic stress release and its subsequent postseismic deformation can be found in the Sumatran-Andaman region, Indonesia (McCloskey et al., 2005; Pollitz et al., 2006; Hughes et al., 2010; Qiu and Chan, 2019) and earthquakes in strike-slip fault systems such as the 1999 Hector Mine earthquake in Southern California, United States (Freed and Lin, 2001) and the 2016 Visso earthquake in Italy (Tung and Masterlark, 2018).

On the 27th of February 2010, the Maule earthquake occurred in Central Chile with a magnitude of Mw 8.8, breaking a 500 km long plate boundary segment along strike (Figure 1 and Vigny et al., 2011; Moreno et al., 2012). About 5 years later, a segment approximately 200 km north of the Maule rupture area ruptured with the Mw 8.3 Illapel earthquake that occurred on September 16th, 2015 (Figure 1). This earthquake broke a presumably mature segment, where the last event occurred in 1943 with a magnitude of Mw 7.9 (Beck et al., 1998). In the Illapel segment (enclosed area in Figure 2), an enhancement of landward displacement was observed following the 2010 Maule event (Figure 2 and Klein et al., 2016; Melnick et al., 2017; Ruiz et al., 2016). Using GPS data, Ruiz et al. (2016) revealed an increase of the annual interseismic velocity in this segment of approximately 15% between the 2010 event and the end of 2014, but they did not measure the locking degree directly or the stress changes due to this increase. Melnick et al. (2017) found an increase of 20% in the apparent locking degree by using GPS observations in the period within the 2010 Maule and the 2015 Illapel events. In the study of Melnick et al. (2017), however, the effect of exponential transient components, e.g., viscous relaxation, were not included, as they only focused on the change in the linear trend component of the observations between the 2010 and 2015 events to calculate the locking changes. By contrast, Klein et al. (2016) showed that a possible mechanism to explain such a pattern may be viscoelastic relaxation in the asthenosphere, as well as in a deep, low-viscosity channel along the fault interface, using a linear Burgers rheology model. Despite these advances in understanding the kinematics of these patterns, it remains unclear how the 3D power-law rheology affects the landward postseismic displacements and how stresses were built-up on the fault interface in the Illapel segment due to the coseismic and postseismic deformation of the 2010 Maule event.

Here, we use the 3D geomechanical-numerical model of Peña et al. (2020) and cumulative GPS observations between the 2010 Maule and 2015 Illapel earthquakes to examine the transient deformation and stress changes in the Illapel segment induced by the coseismic and postseismic deformation of the 2010 Maule earthquake. We first investigate the differences in the resulting surface postseismic displacements due to the choice of model rheology and how this affects the calculation of the Coulomb Failure Stress (CFS) changes. Finally, we evaluate if the CFS changes induced by the co- and postseismic effects of the 2010 Maule earthquake could have triggered the 2015 Illapel event.

We use the 3D geomechanical-numerical model of Peña et al. (2020) to investigate both the deformation and stress transient patterns induced in the Illapel segment due to the 2010 Maule event. The model described in Peña et al. (2020) considers the slab and Moho geometries obtained from Hayes et al. (2012) and Tassara et al. (2006), respectively. The model extends 4000 km in the WE direction, 2000 km in the NS direction, and 400 km in the vertical direction (see Peña et al. (2020) for further details). The postseismic deformation following the 2010 Maule event is simulated by combining afterslip and viscous relaxation jointly through a forward modeling approach. To do so, we employ the best-fit simulation of Peña et al. (2020), which considers afterslip on the fault interface, and power-law rheology with dislocation creep processes in the crust and upper mantle (see elastic and creep parameters in Supplementary Tables 1, 2). In this simulation, most of the afterslip occurs at depths <60 km (Figure 2A), while viscous relaxation occurs in the continental lower crust and upper mantle beneath the volcanic arc (see Figure 8A in Peña et al., 2020). As discussed by Peña et al. (2020), this model was favored not only because of the best fit to cumulative GPS observations and GPS time-series but also because of the correlation between the modeled afterslip and the aftershocks. For comparison, we also consider a simulation with linear rheology, which consists of a whole elastic crust and linear viscosity in the upper mantle. The afterslip distribution from the linear rheology case is displayed in Figure 2B.

In this study, we compute the CFS changes on the megathrust fault interface in the Illapel segment due to the coseismic and postseismic deformation of the 2010 main shock. We calculate the CFS changes employing the linear relationship between shear and normal stresses given by,

where Δτ is the change in shear stresses, μ′ the effective coefficient of friction and Δσ the change in normal stresses (e.g., King et al., 1994; Stein, 1999). Following the approach of Moreno et al. (2018), we also consider a heterogeneous effective coefficient of friction on the fault interface. The study of Moreno et al. (2018) showed that this approach could explain the occurrence of moderate-size earthquakes in southern Chile better. Moreno et al. (2018) estimated μ′ = 0.04 for regions where the locking degree is ≥ 0.8 and μ′ = 0.01 for regions with a locking degree < 0.8. We employ a similar approach, testing for three different locking degree threshold values (0.8, 0.5, and 0.2), where μ′ is set to 0.04 and 0.01 above and below the locking threshold, respectively. These low values of μ′ are supported by a recent study by Dielforder et al. (2020) and by the inversion of stress orientations of moment tensors from subduction earthquakes as well. This suggests a frictional zonation with a higher effective coefficient of friction in highly locked domains as compared to creeping ones (Hardebeck and Loveless, 2018). We tested end-member μ′ values as well (Supplementary Figure 1), but they do not result in any major CFS changes. This agrees with Pollitz et al. (2006), who found similar stress changes by considering a range of μ′ between 0 and 0.8 to investigate the stress changes in the Sumatra-Andaman region. The locking coefficient is taken from Tilmann et al. (2016) who use an inversion of the interseismic GPS observations using an inversion method based on FEM-generated Green’s functions. Their implemented fault geometry and elastic parameters follows the approach of Moreno et al. (2012), which are also consistent with those used in the 3D model of Peña et al. (2020). The resulting numerical problem is resolved with the commercial finite element code Abaqus, version 6.14 and the CFS changes are computed on the fault interface using GeoStress Add-on (Stromeyer et al., 2020) for Tecplot 360 EX, version 2019 R1.

To investigate the differences in deformation of the Illapel segment due to the choices of model rheology, we use the cumulative postseismic GPS displacements from Melnick et al. (2017) between the 2010 Maule and 2015 Illapel earthquakes. In this data set, the effect of aftershocks and seasonal components have been removed (further information can be found in supporting information in Melnick et al., 2017). Following the approach from previous studies (e.g., Klein et al., 2016; Weiss et al., 2019; Peña et al., 2020), we have also removed the interseismic trend measured prior to the 2010 Maule earthquake to thus account only for the postseismic effects following the 2010 Maule event (Figure 2).

The surface displacements are presented in Figure 2 and the coseismic and postseismic CFS changes are presented in Figures 3, 4. Finally, we compare the cumulative CFS changes with the seismicity between the 2010 Maule and 2015 Illapel events, including the Illapel main shock.

Landward deformation patterns following megathrust earthquakes have been observed by modern GPS data (e.g., Heki and Mitsui, 2013; Melnick et al., 2017), but their controlling mechanism is not fully understood. Our results show that the choice of rheology has a large effect on the azimuth of deformation at segments adjacent to the mainshock. Here, an upper mantle with homogenous and linear viscoelastic material properties cannot predict the observed landward pattern. In contrast, a model considering second-order, power-law rheology in the crust and upper mantle can better explain this pattern (Figure 2A). An alternative approach to explain this pattern is the one of Klein et al. (2016), who considered an upper mantle with Burgers rheology, as well as a deep subduction channel with low viscosities (∼10¹⁷ Pa s). On the other hand, Yuzariyadi and Heki (2020) showed that a generic 3D model with linear Maxwell rheology in the upper mantle exhibited trenchward patterns. Because of this discrepancy with the GPS data, Yuzariyadi and Heki (2020) hypothesized that deep slab acceleration may also play a key role, in agreement with the study of Heki and Mitsui (2013). Nevertheless, our findings indicate that the surface displacement patterns from the 3D forward model of Yuzariyadi and Heki (2020) are produced by the simplicity of the rheological structure of their model, which considers a linear viscoelastic mantle, whose viscosity is homogeneously distributed in the whole upper mantle, and a fully elastic crust.

Even though the fit of the postseismic displacement with the forward model that uses non-linear rheology is much better in comparison with linear viscoelastic implementation, it cannot perfectly predict the north component of the GPS signal. However, this has hardly any impact on the estimates of the CFS changes as the stress traction vector is resolved in the rake vector direction of 83° from the Illapel earthquake, which means that there is marginally any contribution from the north component of the displacement to the stress changes.

The CFS changes produced by the coseismic and postseismic deformation of the Maule event are relatively small. Nevertheless, these small values are sufficient to trigger seismicity, as statically shown in other subduction and strike-system settings with a threshold magnitude as small as 0.1 bar (e.g., Stein, 1999, and references therein; Ziv and Rubin, 2000). Moreover, laboratory-derived, rate-and-state friction laws predict that small CFS changes can promote the increase of seismicity (e.g., Stein, 1999, and references therein). Notably, in the southern region R1 (Figure 4A), five events Mw ≥ 5.5 occurred in regions with a change of CFS > 0.1 bar (Figure 4D). Our findings reveal that at the hypocenter of the two larger events (Mw 6.4 and 6.7) the CFS increased by ∼0.63 and 0.14 bar, respectively. Furthermore, our results indicate a potential triggering of the events in region R1 due to the joint contribution of coseismic and postseismic deformation (Figure 4A).

Yet, no strong increase in seismicity rate changes has been reported after the 2010 Maule event in the whole Illapel segment (Supplementary Figure 7 and Heidarzadeh et al., 2016; Poli et al., 2017; Ruiz et al., 2016), but so far only events ≥ Mw 4.0 have been considered. Huang and Meng (2018) used data provided by the Centro Sismológico Nacional de Chile (CSN) to extract repeating earthquakes (Mw ≥ 2.5), employing a matched-filter-enhanced method. They considered ∼4 years before the Illapel earthquake and found an increase of seismic activity in the Illapel segment 100 days before the Illapel main shock but it remains unclear how seismicity evolved in the first ∼2 years following the 2010 Maule event, a period in which the postseismic deformation associated to the Maule event prevails (e.g., Bedford et al., 2016; Klein et al., 2016; Peña et al., 2020). As shown by Toda et al. (2011) for the case of the enhanced seismicity rates following the 2011 Tohoku-Oki earthquake in Japan, an enhanced seismic catalog considering M ≤ 4.0 events may provide key information to ultimately evaluate the triggering of the events in the southern region under the CFS criteria.

On the other hand, at the hypocenter of the Illapel main shock and events further north, the calculated CFS changes (<0.06 bar) are below the assumed threshold CFS of 0.1 bar. Furthermore, if we use the shear stress drop of 30–60 bar of the Illapel earthquake estimated by Yin et al. (2016) with a 3D kinematic model and assuming that this shear stress has been accumulated since the previous Mw 7.9 event on the same segment in 1944 (Beck et al., 1998), the annual shear stress loading rate would be between 0.42-0.84 bar/yr. Changes in the normal stress due to overburden can be assumed to be very small. Consequently, the CFS changes of ∼0.06 bar that is transferred co- and postseismically to the Illapel hypocenter represent only 7-14% of the annual loading rate. Thus, given that the CFS changes are below the threshold of 0.1 bar and the small clock-advance of approximately a few months only, we have two options how to interpret the results. Either there is no triggering and earthquake occurrence is a Poisson process in subduction zones, or another mechanism is acting that is not represented in the model used for the stress transfer calculations.

Additional processes transferring CFS more efficiently may be related to dynamic triggering due to the passage of seismic waves (e.g., Gomberg et al., 2001; Freed, 2005; van der Elst and Brodsky, 2010), or changes in pore pressure due to changes in the permeability in the fault interface (e.g., Brodsky et al., 2003; Elkhoury et al., 2006; Manga et al., 2012). The first one may change the state of the stresses of the medium due to the passage of the seismic waves, both at near and remote distances, but is a very short-term mechanism (e.g., Freed, 2005; van der Elst and Brodsky, 2010). Given the small coseismically CFS values inferred at the hypocenter of very short-term events E1 (0.003 bar) and E2 (0.07 bar, Figure 4A), 9 days and 30 min after the Maule event, respectively (Supplementary Figure 1), dynamic wave processes are likely candidates for the occurrence of these events. On the other hand, pore-fluid processes last longer and produce surface deformation patterns confined close to the rupture zone, i.e., mostly offshore at subduction zones, rather than at further distance (e.g., Masterlark, 2003; Hu et al., 2014; McCormack et al., 2020). Particularly, pore pressure changes would significantly affect the effective normal stresses (e.g., Brodsky et al., 2003; Elkhoury et al., 2006; Hughes et al., 2010; Manga et al., 2012), and would also explain the delay in response due to the slow diffusion of the fluids. In addition, its resulting surface displacement field is generally smaller than that from viscous relaxation, although the CFS changes are larger than those from viscous relaxation only (Masterlark et al., 2001; Hughes et al., 2010). Therefore, our CFS calculations in the Illapel segment during the postseismic period may be underestimated by not considering pore-pressure processes. Its integration into the 3D model would require further constraints from geodetic observations, potentially from InSAR data using its high spatial resolution.

We examine the co- and postseismic deformation and stress change pattern on the Illapel segment induced by the 2010 Maule earthquake using a 3D geomechanical-numerical model with linear and power-law rheology. Our results show that the postseismic horizontal deformation patterns in the Illapel segment are particularly sensitive to the choice of the rheology model. Here, the power-law rheology model can better explain both the landward directions and magnitude of the cumulative GPS displacement vectors. The opposite surface displacement patterns from the choice of rheology also have a major impact on the resulting postseismic CFS changes, but compared to the co-seismic stress transfer, only 15% of the CFS changes are transferred postseismically to the hypocenter location of the 2015 Illapel earthquake. From magnitudes of CFS changes it seems quite unlikely that the Illapel earthquake was triggered by the Maule event - such a small CFS change of <0.06 bar at the hypocenter converts to a clock-advance of approximately 2 months. This finding would not change with a better fit of the modeled north component of the GPS observation. Hence, either no triggering occurred or a local secondary process such a pore-pressure changes, which is not considered here, potentially contributed to a remote triggering mechanism.

We also compared the cumulative CFS changes to the seismicity that occurred between the Maule and Illapel events, using the model with the power-law rheology. Here, particularly in the southern region of the Illapel segment, sixteen seismic events occurred in regions of positive CFS changes above the triggering threshold of 0.1 bar, indicating a potential triggering because of the Maule event. Therefore, our study emphasizes the importance of using power-law rheology models to also investigate surface landward patterns and stress changes induced by megathrust earthquakes in neighboring segments. Future investigation should especially focus on assessing the transfer of stresses from the 2014 Iquique in its southern segment (Loa segment) and from the 2015 Illapel earthquake in its northern segment (Atacama segment), gaps remaining with potential magnitudes exceeding M 8.0.

Publicly available datasets were analyzed in this study. This data can be found here: Melnick D. et al. (2017). “The super-interseismic phase of the megathrust earthquake cycle in Chile.” Geophysical Research Letters 44 (2): 784–791. https://doi.org/10.1002/2016GL071845.

CP: conceptualization, methodology, software, validation, formal analysis, investigation, writing – original draft, writing – review & editing, visualization, supervision, project administration, and funding acquisition. OH: conceptualization, methodology, formal analysis, investigation, resources, writing – original draft, writing – review & editing, supervision, project administration, and funding acquisition. MM: conceptualization, writing – review & editing. DM: conceptualization, formal analysis, investigation, writing – review & editing. OO: conceptualization, writing – review & editing. All authors contributed to the article and approved the submitted version.

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.