# Early Stage and Main Ruptures of the 2015 Mw8.3 Illapel, Chile, Megathrust Earthquake: Kinematic Elliptical Inversions and Dynamic Rupture Simulations

Hideo Aochi, Sergio Ruiz

► To cite this version:

Hideo Aochi, Sergio Ruiz. Early Stage and Main Ruptures of the 2015 Mw8.3 Illapel, Chile, Megathrust Earthquake: Kinematic Elliptical Inversions and Dynamic Rupture Simulations. Journal of Geophysical Research: Solid Earth, 2021, 126 (5), 10.1029/2020JB021207. http://dx.doi.org/10.1029/2020JB021207

HAL Id: hal-03878014

https://hal.science/hal-03878014

Submitted on 29 Nov 2022

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
JGR Solid Earth

RESEARCH ARTICLE

10.1029/2020JB021207

**Key Points:**

* Rapid kinematic slip inversion identifies a main patch in the north preceded by a small one in the south, with a time delay

* Dynamic rupture simulations constrained by kinematic inversion reconstruct the early and main rupture process of the 2015 Illapel earthquake

* The ruptured area and heterogeneous slip patches can be reconstructed from the interseismic coupling map and agree with recorded seismicity

**Supporting Information:**

Supporting Information may be found in the online version of this article.

**Correspondence to:**

H. Aochi,
aochi.hideo@gmail.com

**Citation:**

Aochi, H., & Ruiz, S. (2021). Early stage and main ruptures of the 2015 Mw8.3 Illapel, Chile, megathrust earthquake: Kinematic elliptical inversions and dynamic rupture simulations. *Journal of Geophysical Research: Solid Earth, 126*, e2020JB021207. https://doi.org/10.1029/2020JB021207

Received 22 OCT 2020
Accepted 15 APR 2021

# Early Stage and Main Ruptures of the 2015 Mw8.3 Illapel, Chile, Megathrust Earthquake: Kinematic Elliptical Inversions and Dynamic Rupture Simulations

Hideo Aochi¹,² ID and Sergio Ruiz³ ID

¹Laboratoire de Géologie, Ecole Normale Supérieure, CNRS UMR 8538, PSL Research University, Paris, France, ²Bureau des Recherches Géologiques et Minières, Orléans, France, ³Department of Geophysics, Universidad de Chile, Santiago de Chile, Chile

**Abstract** We apply kinematic and dynamic modeling to the 2015 Mw8.3 Illapel, Chile, earthquake constrained by continuous high-rate GPS and strong motion data. Kinematic inversion by elliptical patches allows us to rapidly outline the ruptured area with different time windows and frequency ranges. The preferred solution indicates that the main large patch in the north is preceded by a small one in the south close to the hypocenter with a time shift no greater than 20 s. The rupture directivity on the main patch suggests that the origin is at depth, away from the initial small hypocentral patch. We then construct dynamic rupture models to be consistent with the geometry obtained from the kinematic inversion. We estimate the fracture energy of the main patch to be ~7.5 MJ/m². The initial rupture of the small hypocentral patch does not trigger the main patch due to the gap in fracture energy, and thus another nucleation is set at depth. This can be regarded as a foreshock-mainshock sequence rather than a direct cascade rupture growth. The ruptured area corresponds to the preexisting zone of large interseismic coupling prior to the Illapel earthquake. The historical seismicity of the previous century suggests a possible reconstruction of the asperity map, assuming that every earthquake represents a characteristic seismogenic patch. Therefore, the construction of dynamic ruptures with geodetic and seismological knowledge is possible and useful not only for reproducing known earthquakes but also for providing a physically constrained model for quantitative seismic hazard study.

**Plain Language Summary** The 2015 Mw8.3 Illapel (Chile) earthquake is considered one of the mega-earthquakes of this century. This study constructs a dynamic rupture model which is consistent with the kinematic inversion and frictional parameter estimates. Our results support that the Illapel earthquake consists of a main large patch preceded by a smaller one by dozens of seconds. The ruptured area is consistent with the interseismic coupling map and the seismicity. These elements suggest that the Illapel earthquake occurred in a comprehensive way, and the shown strategy of simulation is applicable to other earthquakes to improve our ability to construct a constrained mechanical model for quantitative seismic hazard assessment.

## 1. Introduction

The Mw8.3 2015 Illapel (Chile) earthquake had a rupture of ~200-km length and 100-km width from the deepest zone of the plate interface to the trench (Melgar et al., 2016; Ruiz et al., 2016; Tilmann et al., 2016; Ye et al., 2016; among others), as Figure 1 shows the area of our modeling and the related observations. This event was located in north-central Chile, where big earthquakes occur regularly (Ruiz & Madariaga, 2018). Previous events occurred here in 1943, 1880, and 1730; however, these rupture extents were not always the same (Beck et al., 1998; Ruiz & Madariaga, 2018; Tilmann et al., 2016). For example, the 1730 event had a magnitude of ~9.0 (Carvajal et al., 2017; Udias et al., 2012). Modeling the dynamic rupture of such large subduction earthquakes is one of the most important issues for understanding not only the long-term subduction process but also the severe ground motions observed in cities close to rupture zones (Fernandez et al., 2019) and associated tsunami impacts (e.g., Aránguiz et al., 2016; Fuentes et al., 2017; Heidarzadeh et al., 2016; Lay et al., 2016) for quantitative seismic hazard analysis.

The kinematic rupture process of the 2015 Illapel earthquake has been studied by many groups; some use a combination of teleseismic waveform and geodetic data (i.e., Barnhart et al., 2016; Carrasco et al., 2019;
<figure><img src="image_2.png" /><figcaption>**Figure 1.** (a) Modeled area for the 2015 Mw8.3 Illapel, Chile, earthquake. The epicenter is represented by a star. Triangles and squares show strong ground motion (CSN) and continuous GPS stations (GNSS), respectively. Colored grid points represent the fault plane (220 × 140 km, every 10 km) with the vertical position (depth) in color. The local coordinates ($\xi_1, \xi_2$) are defined with respect to the hypocenter position, letting $\xi_1$ and $\xi_2$ point in the updip and along-strike (N185°E) directions, respectively. Small circles show the spatiotemporal position of the high-frequency generation sources (Ruiz et al., 2016). (b) The accelerations of the ground motions in EW component along the longitude at the stations indicated in (a). The figure is modified after Ruiz et al. (2016). The gray shadows indicate the early burst of the energy in the ground motions.</figcaption></figure>

Grandin et al., 2016; Klein et al., 2017; Lange et al., 2016; Meng et al., 2018; Okuwaki et al., 2016; Williamson et al., 2017; Zhang et al., 2016) and the others include tsunami information (i.e., An et al., 2017; Li et al., 2016; Satake & Heidarzadeh, 2017; Williamson et al., 2017). However, the strong motion data of the National Seismological Center and the continuous recording of 1-s sampling of GNSS data (Baez et al., 2018; Leyton et al., 2018a) have not yet been sufficiently explored. As already pointed out by Ruiz et al. (2016) and Okuwaki et al. (2016), the strong motions during the 2015 Illapel earthquake showed multiple bursts of high-frequency energy (Figure 1). A large moment release should have been brought by the wave burst of a large amplitude and a long duration, observed on all the stations. Furthermore, a small, early part of the wave phases is observed particularly at the stations near the epicenter, indicated by gray zone, for the beginning of the records. This early part corresponds to a Mw6.9 subevent according to the estimation of Ruiz et al. (2016). Okuwaki et al. (2016) point out the two distinguishable episodes from the arrivals of *P* and *S* waves in the near station CO03. The spatiotemporal evolution of high-frequency sources from the back projection analyses (also in Figure 1) also indicates a complex rupture trajectory (e.g., An et al., 2017; Meng et al., 2018; Ruiz et al., 2016). The first trajectory progresses from the epicenter to NW direction. On the other hand, a downdip high-frequency radiation patch emerges around 20 s or later, indicating a possible second nucleation patch at depth (An et al., 2017; Meng et al., 2018; Yin et al., 2016). Thus, one can distinguish the first from the second events such as prior and main rupture processes. The multiple bursts in seismograms are similarly observed in the records of the Maule 2010 and Tohoku-Oki 2011 earthquakes, indicating the complex rupture process (Kurahashi & Irikura, 2011; Ruiz et al., 2012). A cascading rupture has been also reported for the 2020 Mw6.8 Elazığ (Turkey) earthquake and simulated dynamically (Gallovič et al., 2020). We now aim to model the 2015 Illapel earthquake in low frequencies within this scope to explore possible dynamic processes for rupture growth: a single rupture growth or multiple-growths (foreshock-mainshock sequence) as interpreted in the seismograms.
In this paper, we first attempt to distinguish the two rupture processes through kinematic inversion. Our approach considers an elliptical patch model, which simplifies the rupture process with very few parame-
ters to rapidly outline the rupture process (rupture area, rupture directivity) as proposed by Vallée and Bou-
chon (2004) and applied for various earthquakes not only kinematically but also dynamically to estimate
the frictional parameters (Di Carli et al., 2010; Peyrat et al., 2010; Ruiz & Madariaga, 2011; Ruiz et al., 2017;
Twardzik et al., 2014; Ulrich & Aochi, 2015). Description by slip patches is consistent with the mechanical
multiscale heterogeneous model of fault (e.g., Ide & Aochi, 2005) and the model of strong motion genera-
tion areas (e.g., Irikura & Miyake, 2011). Previously, the 2011 Tohoku earthquake was dynamically simu-
lated using these concepts from rupture growth (Ide & Aochi, 2013) and seismic wave radiation at different
frequencies (Galvez et al., 2014). It is essential that a particular patch size is predominant for frequencies
of interest. The second objective of this paper is thus to construct a simple, representative dynamic rupture
model. We attempt to estimate the frictional parameters of the dominant patch through dynamic rupture
simulations using the geometry found by the preceding kinematic inversion. Finally, we compare our patch
model with the seismicity and interseismic coupling map to discuss how we can construct a mechanically
understandable model of this megathrust earthquake.

## 2. Data and Model

The September 16, 2015 Illapel earthquake was observed by wide seismological and geometrical networks in Chile and surrounding region. Continuous high-rate GPS observation points are quite uniform from north to south along the strike, while the strong ground motion stations are located mainly in the far North, with few stations being close to the epicenter (Figure 1). We then use mostly cGPS data for their coverage, which is sufficient when using low frequencies up to 0.1 Hz. For the comparison of dynamic simulations, we also use acceleration data for higher frequencies. In accordance with the Centro Sismológico Nacional of the Universidad de Chile (http://www.sismologia.cl/), we set the position of the hypocenter at (W71.741°, S31.637°, 23.3-km depth) and the origin time at 22:54:31 (UTC) for reference. All data and simulations are hereafter aligned at this origin time, noting that a leap time of 17 s is included in the cGPS data and, confirmed by the synchronization between CO06 (acceleration) and PFRJ (cGPS) collocated stations (Baez et al., 2018; Leyton et al., 2018a). The fault plane is set with a strike of N5°E, a dip of 19°, and an assumed rake of 90° for the simplicity of the calculation, which is close to the reported values by different seismo-logical agencies. The fault plane covers an area of 220 km (strike) × 140 km (dip), reaching the trench (Figure 1).

Let us describe the earthquake with elliptical patches, each of which generally requires only eight param-
eters in kinematics (Figure 2a). Two parameters are needed for its position (X<sub>asp</sub>, Y<sub>asp</sub>) and three for its
geometry (R<sub>a</sub>, R<sub>b</sub>, θ). A slip distribution on the patch is assumed of elliptical shape with crack-like slip
distribution with maximum (D) (e.g., Ruiz & Madariaga, 2013; Ruiz et al., 2017; Twardzik et al., 2014; Ul-
rich & Aochi, 2015). The rise time (T<sub>R</sub>) and rupture velocity (V<sub>R</sub>) are set constant on the patch. Compared
to the kinematic descriptions by small subfaults, this parametrization significantly decreases the number
of parameters. We also consider two more parameters, delay time (t<sub>0</sub>) and rupture directivity direction (η),
for a more general imaging of the rupture propagation on a patch, which is independent from the supposed
hypocenter and origin time. Therefore, we simultaneously search 10 parameters for one patch and attempt
to identify up to two patches, sequentially.

We then use the kinematically found patch geometry for the dynamic rupture simulation. In this case, the fault interface is governed by the slip-weakening friction law (Ida, 1972) with peak and residual strengths (τ<sub>p</sub>, τ<sub>r</sub>) and critical slip displacement (D<sub>c</sub>). The fracture energy (G) is defined as G = (τ<sub>p</sub> − τ<sub>r</sub>) × D<sub>c</sub> / 2. In the case where G<sub>c</sub> is scale-dependent (e.g., Ohnaka, 2003), an earthquake rupture is expected to grow from a small nucleation patch to a large patch as a cascade (Ide & Aochi, 2005) (Figure 2b). We will search for a better set of the frictional parameters in the assumed range of the parameters and then discuss possible scenarios for rupture growth.

Both for kinematic and dynamic simulations of this earthquake, the seismic wave propagation is calculated
using the Finite Difference Method (FDM; Aochi & Madariaga, 2003; Aochi et al., 2013) with fourth-or-
der staggered grids. A 1D layered structure (Figure S1) is adopted based on the Crust 1.0 model (Wang
<figure><img src="image_3.png" /><figcaption>**Figure 2.** (a) Definition of model parameters for kinematic source description by patches. See also Table S2. The rupture time is calculated not from the hypocenter but from the rupture initiation position, inclined by an angle of $\eta$, by keeping the same distance to the center of ellipse ($R_{hypo-ellipse}$). Here, the definition of rupture directivity direction, $\eta$, is a difference with respect to the hypocenter-elliptical center direction. (b) Slip-weakening friction defined and schematic illustration of multiscale patch model. $\tau_p$ and $\tau_r$ are peak and residual strengths, respectively, and $\tau_0$ is initial shear stress. The critical slip displacement, $D_c$, is considered proportional to the patch size, $R_n$.</figcaption></figure>

et al., 2003). This is sufficient for comparisons of ground motions in low frequencies, as we principally use frequencies lower than 0.1 Hz. The grid size and time step for this calculation are set as 500 m and 0.02 s, or 200 m and 0.01 s, corresponding to the maximum frequency of 0.43 or 1.07 Hz, respectively.

# 3. Kinematic Inversion by Elliptical Patch Approach

## 3.1. Identification of Main Patch

First, we aim to detect a main elliptical patch from the regional 10 cGPS stations (Figure 1). Both for the synthetic and observation waveforms, we use a digital filter (Saito, 1978). For obtaining 10 unknown independent parameters, we adopt a Genetic Algorithm (GA, e.g., Goldberg, 1989), technically the GA Utility Library (GAUL), following the approach of Ulrich and Aochi (2015). The model parameters are summarized in Table S1. We vary the time window length between 60 and 200 s since we are interested in how quickly the patch (area, rupture propagation, and corresponding magnitude) can be identified. Figure 3 shows an example of the inversion process for a time window of 180 s and a frequency band between 0.01 and 0.05 Hz. The GA consists of generation and population to stochastically select the individuals for a better fit. The misfit ($\lambda$) is calculated using the squared difference between the observed and synthetic waveforms for the three components

$$ \lambda = \sum_{\text{component}} \sum_{t} (Y_{\text{syn}}(t) - Y_{\text{obs}}(t))^2 \quad (1) $$

After several attempts, we decide not to give a weight in Equation 1, as the inversion focuses on the macroscopic feature of the rupture process. We run the GA algorithm until the 21st generation. The main findings are that the patch center is found in the North around (72°W, 31°S), the ruptured area reaches the trench, and the rupture propagation is delayed with respect to the given origin time, starting from the deepest part.

We run 10 random inversion processes at each time window length. We show a typical geometrical solution at every time window length illustrating how they differ (Figure 4a) and the mean and standard deviation
<figure><img src="image_4.png" /><figcaption>**Figure 3.** Inversion process for a time window of 180 s. (a) The convergence of misfit and 10 inverted parameters. (b) The slip distribution and rupture time from the obtained solution. Star is the epicenter location. Cross represents the rupture initiation position defined by $t_0$ and $\eta$. The triangles show the positions of the cGPS stations. For reference, the Centroid-Moment-Tensor of Global CMT is also shown. (c) The comparison between the observation (black) and the final synthetics (red). The intermediate solutions at every generation are shown by gray lines.</figcaption></figure>

for the obtained magnitude, patch size, and fault slip (Figure 4b). At each time window length, the solutions are numerically converged, as demonstrated in Figure 3; however, we need further interpretations to compare the solutions of different time windows. At shorter time windows, only a portion of the patch is ruptured, or the seismic wave radiation has not yet arrived at the far stations (rupture time on the patch ranges between 20 and 100 s in Figure 3). Each inversion reflects the characteristics of the rupture process according to the arrived phase at a given time window, and we need to determine if the solutions are stable or in transition. The solution should be more reliable when using a longer time window, and the estimated magnitude should become stable after 140 s (Figure 4b). The estimated magnitude reaches around Mw8 for a time window of 150 s, and the solutions become stable (Mw7.99, ruptured area $1.3 \times 10^4$ km$^2$, and maximum slip 4.7 m averaged for time windows of 150–200 s). The underestimation of magnitude is principally due to the simplification of the parametrization (geometry, slip shape) as well as the limited frequency range of the analysis.

Let us use a low-pass filter up to 0.05 Hz to keep the permanent ground displacement. Figure 5 shows an example of the inversion with a time window of 150 s, and Figure 6 summarizes the solutions for different time window lengths. Because of a wider frequency range, particularly with the permanent displacement, the solutions are well-constrained for quickly overviewing the rupture process, without being influenced by the signals of higher frequencies. The solutions become stable after 100 s (magnitude 8.13, ruptured area $1.4 \times 10^4$ km$^2$, and maximum slip 6.7 m averaged over time windows of 100–180 s. This meets the necessary criteria for showing how early the approximate finite source parameters and magnitude can be estimated to
<figure><img src="image_5.png" /><figcaption>**Figure 4.** Solutions for different time window lengths for the frequency band between 0.01 and 0.05 Hz. (a) The typical geometry of the patches obtained for different time window lengths between 60 and 200 s from the preset origin time $t = 0$. (b) The moment magnitude, the geometrical average of elliptical semiaxes, and the fault slip at the center of the ellipse as a function of the time window length. The values are averaged after 10 runs at each time window, and error bars indicate standard deviation.</figcaption></figure>

provide the first alert of subduction earthquakes (Leyton et al., 2018b). Besides of Equation 1, we calculate a variance reduction defined as

$$VR = 1 - \frac{\sum (Y_{\text{syn}}(t) - Y_{\text{obs}}(t))^2}{\sum Y_{\text{obs}}(t)^2} \quad (2)$$

and it is confirmed that the use of low-pass filter leads to a stable result (VR close to 0.9) (see Figure S2).

From these tests, the finite source parameters of the main patch can be found without waiting for the passage of all the waves. In other finite source inversions by subfaults (e.g., An et al., 2017; Melgar et al., 2016), the peak of source-time function is found for 50–60 s. Thus, a time window of a minimum of 100 s is sufficient to detect radiating seismic waves at nearby stations. The identified main patch fits the brief feature of waveforms and permanent displacements. The fact that our large patch is found near the Centroid-Moment-Tensor solution guarantees the reliability of the solution and also infers a rupture propagation toward this position from the epicenter area (e.g., Zahradník et al., 2008). However, we also find that the wave phases at higher frequencies are present at the beginning of waveforms (<100 s at the near stations CNBA, SLMC, LVIL in Figure 3). Next, we investigate this early waveform phase.

## 3.2. A Patch for Early Stage of Waveforms

The visible early phase of waveforms can be treated directly, but we subtract the first inversion result from the full waveforms before the second inversion. We do not simultaneously inverse two patches because this works well only if two patches have similar weights (similar sizes) on the observation (Di Carli et al., 2010; Peyrat et al., 2010; Twardzik et al., 2014; Ulrich & Aochi, 2015). The deconvolution process is similar to the analysis of multiple source-time functions or multiple moment tensors (e.g., Kikuchi & Kanamori, 1991; Tsai et al., 2005; Zahradník & Sokos, 2014). Our approach provides finite source parameters, in particular the finiteness of the rupture area. To capture the early phase of seismic waves, we use a higher frequency band (0.03–0.1 Hz in velocity) from the five closest stations and limit the source area by 45 × 45 km around the hypocenter, for which the Greens functions are calculated every 3 km. We invert nine parameters: ellipse's location ($X_{asp}$, $Y_{asp}$), semiaxes ($R_a$, $R_b$), orientation $\theta$, maximum fault slip $D$, rise time $T_R$, and rupture
<figure><img src="image_6.png" /><figcaption>**Figure 5.** Inversion process for a time window of 150 s and low-pass filtered data at 0.05 Hz. (a) The convergence of misfit and 10 parameters. (b) The slip distribution and rupture time from the obtained solution. Star and cross represent the epicenter location and the rupture initiation position defined by $t_0$ and $\eta$. The Global CMT solution is illustrated for reference. (c) The comparison between the observation (black) and the final synthetics (red). The intermediate solutions at every generation are shown by gray lines. The vertical dotted lines represent 60, 120, and 180 s.</figcaption></figure>

velocity $V_R$. The time lag of the rupture is initially set to zero ($t_0 = 0$), and a mutation process of GA allows for possible variations.

For the deconvolution, we use two models on the main patch: A solution obtained from bandpass filter analysis (time window of 180 s, shown in Figure 3), and another from low-pass filter analysis (time window of 150 s, shown in Figure 5). Figures 7a and 7b show an example of inversion for a time window of 60 s using the former model for deconvolution. The amplitude of the waveforms is the largest for the CNBA station during this early period, signifying the existence of a patch close by. Stations further south (ZAPA, SLMC, C23O) have very low amplitudes and are not as greatly affected by the early rupture process. It is thus understandable why the obtained patch is located in north close to the CNBA station.

We vary the time window lengths for both models of the main patch (Figures 7c and 7d). The patch geometry and corresponding magnitudes are summarized as a function of time window length. When a large phase appears at 30–50 s at the CNBA station, the inversions including this period provide a stable solution, signifying that there is a stationary stage before the growth of the main rupture. In particular, the patch geometry can be identified for the solutions of the time window of 50 and 60 s, regardless of the two
<figure><img src="image_7.png" /><figcaption>**Figure 6.** Solutions for different time window lengths for low-passed data at 0.05 Hz. (a) The typical geometry of the patches obtained for different time window lengths between 60 and 200 s. (b) The moment magnitude, the geometrical average of elliptical semiaxes, and the fault slip at the center of the ellipse as a function of time windows. The values are averaged over 10 runs, and the error bars represent standard deviation.</figcaption></figure>

deconvolutions. The size of this patch remains relatively small—<20 km. Thus, the early process of the rupture can be interpreted to be around Mw6.8, which is close to the Mw6.9 estimated directly from strong motion by Ruiz et al. (2016). This finding is indeed confirmed in the source-time function (Figure 7e). We affirm that the source-time function of the main patch starts at ~20 s, while the small patch starts before.

# 4. Dynamic Simulations

## 4.1. Model Setting

We have learned from the previous kinematic inversions that there are at least two patches playing significant roles: an initial small patch and a main patch controlling the earthquake size. Here, we attempt to reconstruct dynamic rupture models to glean characteristics of the patch distribution for this earthquake. We use multiple circular patches, whose sizes are discretized by a factor of 2 for simplicity, as schematically shown in Figure 2b (Ide & Aochi, 2005). This simple discretization (Ide & Aochi, 2013) allows us to discuss the correlation between the coseismic rupture pattern, seismicity and scaling law.

The background (fault portion outside of the patches) is assumed to slip freely without any stress changes after reaching the peak strength. Using the previous kinematic inversions, a dimension of the main patch ($\sqrt{R_a R_b}$) is ~60–70 km (Figures 3–6). We adopt a radius of 60 km, and then discretize other patch sizes by a factor of 2, namely a radius of 30, 15, and 7.5 km, named Rank 1–4, respectively. We set the main patch (Rank 1) with an elliptical shape to be consistent with the kinematic model as shown in Figure 8a. The principal model parameters are summarized in Table S2. The scaling relation in fracture energy is summarized in Aochi and Twardzik (2020). The small and large patches are slightly superposed, but it is essential that a small patch of Rank 3 does not directly trigger a large patch of Rank 1 due to an important gap in fracture energy. Therefore, it is necessary to introduce a transition process: A medium-size patch (Rank 2) connects the two patches, allowing a cascade rupture growth (Figure 8b), or a second nucleation is prepared at depth at 20 s and grows differently (Figure 8c).

We apply the boundary integral equation method (BIEM; Aochi et al., 2000) in a 3D infinite, homogeneous medium with a mirror source approximation of the free surface effect (Ide & Aochi, 2013). This simplified simulation is suitable for running many simulations. The focus of this paper is not the dynamics of the surface rupture at the trench but the seismic wave radiation from the slip patch at depth. The grid size is 1 km, and the time step is 0.081 s, after a verification test (Figures S3–S7), we also use a coarser grid of 1.5 km and
<figure><img src="image_8.png" /><figcaption>**Figure 7.** Second inversion results for the early process of the Illapel earthquake after the deconvolution of the source models found in the first inversion (a solution from a bandpass filter of 0.01–0.05 Hz with a time window of 180 s and from a low-pass filter up to 0.05 Hz with 150-s window length. (a) The convergence of each parameter for a time window of 60 s after the deconvolution of the latter model. (b) Corresponding fitting of the waveforms between 0.03 and 0.1 Hz for the selected four stations for the parameters shown in (a). (c) The geometry of the patches obtained for different time window lengths between 20 and 70 s. The solid and broken lines correspond to the two deconvolutions before this inversion. Star shows the epicenter, and Global CMT solution is illustrated for reference. (d) The moment magnitude, the geometrical average of elliptical semiaxes, and the fault slip at the center of the ellipse as a function of time windows between 20 and 70 s. The values are averaged over 10 runs for each deconvolution. (e) The source-time function (seismic moment release rate). Two broken lines indicate the two models used for the deconvolution. The solid lines (almost overlapped) show the four inversion results for time windows of 50 and 60 s and both deconvolutions.</figcaption></figure>

time step of 0.121 s to accelerate the parametric study. In each example, a tiny patch (Rank 4) is added at each nucleation point to reassure the rupture growth. Dynamic rupture starts with an instantaneous stress drop on a circular crack of a 6-km radius at a given time with no quasi-static nucleation process.

## 4.2. Synthetic Inversion Test

For the earthquake dynamics, not only the heterogeneity but also the nucleation point is important, as clearly pointed out both kinematically and dynamically by Mai et al. (2005) and Ide and Aochi (2005), e.g. For the patch models presented in Figure 8, the large patch is not simply triggered by the small patch (Model 1). A transition is necessary for the rupture to grow large enough to trigger the large patch. Models 2 and 3 appear to have different timing and directivity of rupture propagation. Here, we carry out a synthetic inversion test on the Models 2 and 3 to show that the kinematic inversion by the patch method
<figure><img src="image_9.png" /><figcaption>**Figure 8.** Three heterogeneous dynamic source models displayed in terms of critical slip displacement ($D_c$) for rupture simulations of the 2015 Illapel earthquake. Patches of four different sizes are distributed. (a) A large patch ($R = 60$ km) and a small patch ($R = 15$ km) are distributed according to the kinematic inversions, with the first nucleation point A. These are the same for all the other models. A tiny patch ($R = 7.5$ km) is given for each nucleation point to secure the rupture growth. (b) An additional medium-size patch ($R = 30$ km) is added. (c) The second nucleation point B is positioned at the edge of the main patch, separately from the small patch. The simulated rupture times and final slip are shown in the middle and at the bottom.</figcaption></figure>

is reproducible. We follow the same procedure as in the previous chapter: namely, the first target is the main patch using 10 cGPS stations with a longer time window (60–180 s) and a lower frequency band (0.01–0.05 Hz), and the second is the early phase of the rupture using the four closest cGPS, with a short time window (20–70 s) and a higher frequency range (0.03–0.1 Hz). The Greens' functions are the same for the synthetic ground motions and inversions.

Figure 9 summarizes the solutions obtained at each step. The main, large patch is found in the North, and a small one is in the South. Only a limited number of stations are shown for the purpose of differentiating the characteristics of the two models, which are clearer in the near-field stations (LVIL in Figure 9b and CNBA in Figure 9c). The geometry of the patches is similar for both models so that it is difficult to distinguish them only from the spatial distribution. The small patch for Model 2 is found slightly further north compared to Model 3. This may indicate that the rupture process of Model 2 is a continuous process leading to the large patch. The shape of the large patch varies with the time windows, and seismic moment is underestimated for Model 3. This is principally because the rupture process is more complex than a simple, uniform kinematic patch. There is a clear difference in the source-time functions: In Model 2, a small patch is a part of the simple bell-shaped source-time function, while in Model 3, the contribution from the two patches can be identified, and there is an important acceleration in the main rupture. Therefore, the small patch in Model 3 can be different from the rupture process of the main patch. As long as the two patches are identifiable in time, they are considered to be mechanically different, such as multiple shocks.

This synthetic test also displays an interesting aspect of the inversion process. For the same time window of 60 s, two different inversion steps provide patch solutions of completely different scales. This is mainly because the stations and phases are preselected for the early-stage inversion. The large-scale inversion is likely
<figure><img src="image_10.png" /><figcaption>**Figure 9.** Synthetic inversion results for target models corresponding to Models 2 and 3 presented in Figure 8. (a) Inferred patch geometry as a function of time window lengths. The broken and solid lines represent two inversion processes of the large scale and the early stage equivalent to the previous analyses. (b) Comparison of ground motions at the selected stations for the large-scale inversion (displacement between 0.01 and 0.05 Hz, time window of 150 s). The gray lines show the intermediate solutions during the inversion. (c) Comparison of ground motions at CNBA (velocity in 0.03–0.1 Hz, time window of 50 s). (d) Comparison of source-time function (rate of seismic moment release $M_0$). All the solutions for different time windows are shown.</figcaption></figure>

to give a large patch solution suitable for low frequencies, even if most part of the patch does not yet contribute to the radiation of seismic waves. Thus, it is important to select the data coverage, frequency bands, and acceptable parameter ranges so as to satisfy the mechanical causality, in particular, if the wave radiation from the obtained ruptured area is included in the waveforms of the analyzed time window length.

### 4.3. Parametric Study on Dynamic Rupture Model

Now, we attempt to find a better parameter set of dynamic rupture simulations. The rupture stability can be discussed in terms of the nondimensional parameter $κ$, introduced by Madariaga and Olsen (2000)

$$ \kappa = \frac{\Delta \tau^2}{\mu \Delta \tau_b} \frac{L}{D_c} \qquad (3) $$

where stress drop is defined as $\Delta \tau = \tau_0 - \tau_r$, with initial stress $\tau_0$ and breakdown strength drop $\Delta \tau_b = \tau_p - \tau_r$ (see also Figure 2b). L represents the characteristic scale of the rupture process, such as asperity size or fault
<figure><img src="image_11.png" /><figcaption>**Figure 10.** (a and b) Two examples of model parameter distribution ($D_c$) and two nucleation points. (c) Misfit as a function of parameter $\kappa$ and initial stress level $\tau_0$. In the gray zone ($\kappa < 0.741$), the large patch does not rupture.</figcaption></figure>

width. For the largest patch of our interest, we test $\Delta\tau = (5, 4, 3.3, 3, 2.5, 1.7)$ MPa, $\Delta\tau_b = (10, 8, 5)$ MPa, and $D_c = (6, 5, 4, 3, 2, 1)$ m. For estimating the parameter $\kappa$, we take $L = 60$ km, which is the diameter of Rank 2 patch, namely, the rupture size when entering the large patch (Rank 1). The misfit is calculated in the same way as in Section 3.1. The first nucleation point is always at the epicenter position. The second nucleation point is set as Model 3 of Figure 9 (Figure 10a) and also at a deeper position (Figure 10b).

Figure 10c summarizes the parametric study as a function of the parameter $\kappa$. For a small $\kappa$ (gray zone), the large patch does not rupture. During the iterations, we also find that a sufficiently small $\kappa$ ($\kappa \sim 1$) leads to a better waveform fit. Thus, we can limit our parameter search without running all the 216 possible parameter sets. The best solution we find is shown in Figure 11a, with $\Delta\tau = 2.5$ MPa, $\Delta\tau_b = 5$ MPa, and $D_c = 3$ m, corresponding to a value of $\kappa = 0.772$, which indicates a mode of subshear rupture (Olsen & Madariaga, 2000; Ruiz & Madariaga, 2013). The order of mean stress drop of 2.5 MPa in our simulations is roughly consistent with the other estimations from finite source models by Yin et al. (2016), regardless of their heterogeneity. The nearest solution (Figure 11b) has the same frictional parameters but a deeper nucleation point. The rupture nucleation time of the second rupture is 20.2 s for the best solution, noting that the smallest hypo- central patch starts at $t = 0$. Such a delay is consistent with the previous results of the kinematic inversions and sufficient for identification using the synthetic inversion test (Figure 9). Finally, we estimate that the fracture energy from the best estimation is $G_c = (\Delta\tau_b D_c)/2 = 7.5$ MJ/m$^2$. It is known that $\Delta\tau_b$ and $D_c$ gen- erally have a bias (e.g., Guatteri & Spudich, 2000), while this might be further distinguished by near-field data (Fukuyama & Mikumo, 2007; Weng & Yang, 2018). Our current approach smoothens any possible heterogeneity inside the patch and is limited in low frequencies. We thus focus principally on the integrated quantity $G_c$. The total fracture energy integrated over the surface is $8.5 \times 10^{16}$ J, for a patch of equivalent radius of 60 km. This estimation consistently meets the dynamic rupture scaling relation compiled by Aochi and Twardzik (2020).

**5. Discussion**

By varying the time window length in the kinematic inversion by the patch method allows us to characterize two patches of different sizes and quickly characterize the finite source parameters without waiting for a chain of seismograms. However, patch description in a short time window may not be perfectly constrained, as it partially contributes to the waveforms by the current time. The need to analyze the obtained solutions remains. It is challenging to characterize the real-time rupture process without using entire seismograms (e.g., Sánchez-Reyes et al., 2018; Uchide & Ide, 2007).
<figure><img src="image_12.png" /><figcaption>**Figure 11.** (a and b) The rupture time and final slip distribution from dynamic rupture simulations for the initial condition of Figures 10a and 10b, respectively. The model (a) minimizes the misfit better than model (b). (c) The comparison of ground motions. The observation is shown by the broken black line. The waveforms are filtered between 0.01 and 0.05 Hz. The two simulations are shifted by 20.2 and 13.8 s, respectively, to minimize the misfit.</figcaption></figure>

The characteristics of slip patches we obtained are consistent with other studies of this earthquake in size, geometry, and timing of rupture process. The irregular rupture process of the Illapel earthquake appears in the change of rupture directivity along the dip (e.g., An et al., 2017; Meng et al., 2018; Ruiz et al., 2016; Tilmann et al., 2016) and the evolution of high-frequency generation areas (Ruiz et al., 2016; Yin et al., 2016) as reviewed in Section 1. Our inversion process is particular in the points that the main patch does not depend on the hypocenter location in its position and rupture directivity, and that fault heterogeneity is simplified enough by definition of the model. This is why our inversion emphasizes the two-step rupture process comparing with the continuous process of the other inversion results.

The difference between the kinematic and dynamic interpretations is that, as long as the causality of rupture time is satisfied, the model description for the former is more flexible, while the latter needs more
<figure><img src="image_13.png" /><figcaption>**Figure 12.** (a) Slip distribution (every 2 m) from the preferred dynamic simulation (Figure 11a). The maximum slip is 7.99 m, and the magnitude is Mw8.16. The earthquakes (stars) are based on the bulletin of International Seismological Center (ISC) for the period 1904–2014 (De Giacomo et al., 2018; ISC, 2018) and are colored according to the following magnitude bins $M_w > 7.7, 7.7 \ge M_w > 7.1, 7.1 \ge M_w > 6.5$. The remaining events are indicated by open stars ($6.5 \ge M_w > 6.0$). Note that the Mw8.1 event occurred on April 06, 1943 and that no earthquakes are reported for the second magnitude bin around the Illapel earthquake. The background represents the coupling obtained by Métois et al. (2016). The color scale is reversed from the original figure to emphasize the strong coupling area with a large coseismic slip. (b) The attributed patch distribution corresponds to each earthquake by color. The shaded area represents area with coupling coefficient smaller than 0.7. The dynamic rupture model obtained in the previous section is illustrated using red lines. (c) Adopted dynamic rupture model for the final simulation in color scale. The main and neighboring patches are shifted to fit the previously obtained patch position. Points A and B denote the first and second nucleation points in the simulation.</figcaption></figure>

constraints on energy balance at every moment. This is why a second nucleation is needed differently from the first nucleation in the dynamic rupture simulations. The dynamic parameters of this mega-earthquake are similar to those of the 2011 Mw9.0 Tohoku (Japan) earthquake (Ide & Aochi, 2013) and the 2016 Mw7.8 Kaikōura earthquake (Ulrich et al., 2019). More recently, Gallovič et al. (2020) obtained for the 2020 M6.8 Elazığ (Turkey) earthquake $D_c \sim 1$ m and a total fracture energy of $3.3 \times 10^{15}$ J. This is close to our scaling assumption (e.g., Aochi & Twardzik, 2020). Yao and Yang (2020) estimated a smaller $D_c$ of 0.25 m and a fracture energy of 0.35–0.7 MJ/m² for the 2012 Mw7.6 Nicoya earthquake. Although the ruptured surface is not precised, we can guess a total fracture energy of $0.35–0.7 \times 10^{16}$ J for the maximum surface of 100 × 100 km. This is relatively small for this size of earthquake, probably because the ruptured area is very localized or this is a very shallow event inferring a very weak fault strength.

In Figure 12, we demonstrate a possible relation between our dynamic rupture model, the regional seismicity and interseismic coupling. We use the catalog from the International Seismological Centre, which provides information from the beginning of the 20th century (Version 5.1; ISC, 2018) for magnitudes larger than 6.5. The localization at the early period may not be precise (Ruiz & Madariaga, 2018), but it allows us to discuss the tendency. The maximum event (Mw8.1) on April 06, 1943 (Beck et al., 1988) corresponds well with the 2015 Illapel earthquake’s location, and they possibly share the same seismogenic elliptical patch. The area has not experienced earthquakes of magnitudes near 7.5, and the second main event is
the Punitaqui earthquake, Mw7.1 on October 15, 1997, an intraplate intermediate depth event of a focal depth of 57.9 km (Lemoine et al., 2001; Pardo et al., 2002). The third group ($7.1 \ge M_w > 6.5$) has several earthquakes, distributed more in the North. The smaller earthquakes occur uniformly along strike, but it seems that they are distributed in the surrounding of the Illapel earthquake as opposed to the center of the ruptured area. This smaller seismicity could define the border of the Illapel earthquake rupture (Poli et al., 2017). We also add the interseismic coupling coefficient of Métois et al. (2016), based on the compiled GPS velocity field prior to the 2015 event. As pointed out in Métois et al. (2016), the ruptured area of the 2015 Illapel earthquake matches well with a place where the seismic coupling is strong and is stopped by low seismic coupling area in the north (~30.5°S). Toward the south (33°S), there have been other earthquakes of magnitudes larger than 7.2 (1971 La Ligua Mw7.8, 1906 Valparaiso Mw8.2). In many cases, the Illapel rupture area is independent from the southern part, except for the huge event of 1730.

Ide and Aochi (2013) propose that the plate interface is characterized by heterogeneous multiscale, intrinsic heterogeneity, which is mapped by patches according to the known seismicity. In Figure 12a, the earthquake locations are colored with different scales, which correspond to different ranks of patch size used in dynamic rupture modeling. We regard the three colored earthquakes (although there is no second group in the area of interest) as having the characteristic circular patch radius of 60, 30, and 15 km (Figure 12b). We neglect the March 28, 1965 Mw7.40 intraplate intermediate depth earthquake (focal depth is 70 km) (Malgrange et al., 1981), since our interest here is limited only on the plate interface. Although the locations of the historical events may be imprecise, we can discuss how the patches identified from the Illapel earthquake correspond to the historical seismicity, as previously discussed for the 2011 Tohoku earthquake (Ide & Aochi, 2013). In other words, we demonstrate here how to construct a dynamic model according to the seismicity. The largest Rank 1 patch of the Illapel earthquake should correspond to the 1943 earthquake (from panels (b) and (c) in Figure 12). Smaller patches are less certain in their positions, but there exists a small cluster of earthquakes with magnitudes smaller than 6.5 (Rank 3) around the epicenter of the Illapel earthquake. We then consider that one or some of them correspond to the initial rupture process of the Illapel earthquake (point A). There is no known earthquake precisely around the second nucleation point, but many earthquakes are mapped around the deeper edge of the Rank 1 patch. Thus, we consider that one of such historical earthquakes initiates the second rupture (point B). The patch distribution of panel (b) can be accordingly adjusted to the positions identified from the Illapel earthquake (Figure 12c). As the point B is slightly apart from the main patch, the nucleation time is set to 16 s.

Figure 13 shows our demonstrative rupture model and compares the ground motions both at cGPS and acceleration stations. We adopt initial stress on Rank 2 patch as $\tau_0 = 2.8$ MPa after testing the values between 2.5 and 4.0 MPa. Although the model remains simple and is only calibrated for low frequencies throughout this study, we intend to show ground motions for higher frequencies at strong motion stations. Our modeling generates some phases seen in the observations related to the change of rupture directivity, rupture timing, and ruptured area consistent with various patches. The numerical simulations reveal the generation of a significant wave when the rupture reaches the trench (the main phase appearing in the EW component of CO06), and this spreads to the other stations. Thus, such dynamic rupture modeling demonstrates its applicability to the broadband ground motion simulations by further calibrating missing details. Therefore, this final demonstration shows that the construction of dynamic ruptures with geodetic and seismological knowledge would be possible with sufficient data sets and useful not only for reproducing known earthquakes but also for providing a physically constrained model for quantitative seismic hazard study.

In this study, we do not explore the mechanical question of how the second nucleation starts. Our focus is rather on the question of how the first rupture ends and the second one is in different dynamic processes. The earthquake initiation is a complex process (e.g., Ellsworth & Beroza, 1995) and always under debate (e.g., Bouchon et al., 2021; Ellsworth & Bulut, 2018). This is technically difficult when different processes occur in a short time span due to the superposition of seismic waves. Delayed rupture process has been reported in different earthquakes, from a delay of about 3 s and a few kilometers away (Abercrombie & Mori, 1994 for the 1992 Landers earthquake) to 70 s and 200 km away (Henry et al., 2000 for the 1998 Balleny, Antarctic, earthquake), e.g. However long-way triggering process is sometimes questioned, as Hjörleisdóttir et al. (2009) propose a continuous rupture process between two subevents for the 1998 Balleny event. Such difference in interpretation is still the case for the 2015 Illapel earthquake. Meng et al. (2018) propose
<figure><img src="image_14.png" /><figcaption>**Figure 13.** (a) Dynamic rupture simulation in terms of rupture time, final slip, and maximum slip rate for the model shown in Figure 12c. (b) Comparison of the ground motions for GPS and strong motion (SM) stations for different frequency ranges. The observations are in black above, and the synthetics are in color at the bottom. The strong ground motion data are integrated once and twice, respectively. The horizontal components of CO06 are exchanged with respect to the originally provided data. The vertical component of CO06 contains an abnormal jump in acceleration, therefore, the integrated waveforms are not reliable for quantitative comparison (not shown).</figcaption></figure>

a bifurcation and detour of the rupture propagation. It can be kinematically possible; however, it is difficult to realize it dynamically. Such detoured rupture should not be visible in wave radiation and final slip, but large enough to continue propagating spontaneously. We have attempted a rupture model connecting two nucleation points; however, this generates an unexpected rupture directivity effect in the NE direction and/or a large slipped area between them. Thus, we have adopted the scenario in which the second nucleation is different from the first. Our numerical simulations indicate that dynamic shear stress change at the second nucleation point is ~0.2 MPa only at 16 s, when the rupture starts. This is still small compared to the fault strength and stress drop, but it is nevertheless possible to dynamically trigger another earthquake (e.g., Voisin, 2002).
## 6. Conclusion

We model the 2015 Mw8.3 Illapel, Chile, earthquake through kinematic and dynamic modeling. First, we perform kinematic inversions by patches using continuous GPS data. We obtain the main patch, primarily responsible for the magnitude of this earthquake. We vary the time window used for inversion. After 90 s from the origin time, the inversion process can outline the finite source parameters of this earthquake, and the ruptured area, which is found in the North, to fit the main phase of the ground motions. The early part of the ground motions (for around 60 s) suggests the existence of a smaller patch; however, the small patch may not have triggered the large main patch. The rupture directivity on the main patch is not from the small patch but from another point located in a different place along plate interface. This two-step rupture process is consistent with the near-field ground motions.

We then construct dynamic rupture models to be consistent with the kinematic inversion. The fracture energy of the main patch is estimated as 7.5 MJ/m². The first rupture initiation, starting at the origin time (t = 0), does not trigger the main patch in the North. Otherwise, the rupture duration is too short and generates a strong rupture directivity to the north. We then set the second rupture initiation at a deeper position on the main patch, and we regard this as a foreshock-mainshock sequence with an interval of 16–20 s. The ruptured area corresponds to the high coupling zone of the interseismic coupling map prior to the Illapel earthquake. The seismicity over a century suggests a possible reconstruction of spatial patch distribution (seismogenic asperities), assuming each patch is intrinsic while leading to characteristic earthquake. This representative dynamic model, accounting for the seismologic and geodetic knowledge, is important not only for reconstructing known earthquakes but also for presenting a physical constrained model for quantitative seismic hazard study.

## Data Availability Statement

The seismicity catalog is from International Seismological Centre (http://www.isc.ac.uk/ and https://doi.org/10.31905/d808b825). The interseismic coupling data (Métois et al., 2016) can be obtained via http://perso.univ-lyon1.fr/marianne.metois/. The used seismological data of Centro Sismológico Nacional of Universidad de Chile can be obtained via http://www.sismologia.cl/. The used continuous GPS data can be obtained via http://gps.csn.uchile.cl/ and the strong motion data from http://evtdb.csn.uchile.cl/. The Global CMT solution (Dziewonski et al., 1981; Ekström et al., 2012) is obtained via https://www.globalcmt.org/. The forward modeling codes of BIEM (Aochi et al., 2000) and FDM (Aochi & Madariaga, 2003) are available at http://doi.org/10.5281/zenodo.1472238 and http://doi.org/10.5281/zenodo.3763864. GAUL is an open source programming library, accessible at http://gaul.sourceforge.net/index.php. Seismic Analysis Code (Goldstein & Snoke, 2005) and Genetic Mapping Tools (https://www Villanueva.org/generic-mapping-tools) are used for signal processing and illustrations (last access October 24, 2019).

## References

* Abercrombie, R., & Mori, J. (1994). Local observations of the onset of a large earthquake: 28 June 1992 Landers, California. *Bulletin of the Seismological Society of America*, 84, 725–734.

* An, C., Yue, H., Sun, J., Meng, L., & Báez, J. C. (2017). The 2015 Mw8.3 Illapel, Chile, Earthquake: Direction-reversed along-dip rupture with localized water reverberation. *Bulletin of the Seismological Society of America*, 107, 2416–2426. https://doi.org/10.1785/0120160393

* Aochi, H., Fukuyama, E., & Matsu'ura, M. (2000). Spontaneous rupture propagation on a non-planar fault in 3-D elastic medium. *Pure and Applied Geophysics*, 157, 2003–2027. https://doi.org/10.1007/PL00001072

* Aochi, H., & Madariaga, R. (2003). The 1999 Izmit, Turkey, earthquake: Nonplanar fault structure, dynamic rupture process, and strong ground motion. *Bulletin of the Seismological Society of America*, 93, 1249–1266. https://doi.org/10.1785/0120020167

* Aochi, H., & Twardzik, C. (2020). Imaging of seismogenic asperities of the 2016 ml 6.0 Amatrice, central Italy, earthquake through dynamic rupture simulations. *Pure and Applied Geophysics*, 177, 1931–1946. https://doi.org/10.1007/s00024-019-02199-z

* Aochi, H., Ulrich, T., Duccellier, A., Dupros, F., & Michea, D. (2013). Finite difference simulations of seismic wave propagation for understanding earthquake physics and predicting ground motions: Advances and challenges. *Journal of Physics: Conference Series*, 454, 012010. https://doi.org/10.1088/1742-6596/454/1/012010

* Aránguiz, R., González, G., González, J., Catalán, P. A., Cienfuegos, R., Yagi, Y., et al. (2016). The 16 September 2015 chile tsunami from the post-tsunami survey and numerical modeling perspectives. *Pure and Applied Geophysics*, 173, 333–348. https://doi.org/10.1007/s00024-015-1225-4

* Báez, J. C., Leyton, F., Troncoso, C., del Campo, F., Bevis, M., Vigny, C., et al. (2018). The Chilean GNSS network: Current status and progress toward early warning applications. *Seismological Research Letters*, 89(4), 1546–1554. https://doi.org/10.1785/0220180011

## Acknowledgments

We thank all our colleagues at ENS Paris and the University of Chile for fruitful discussions and useful information on the Chilean subduction. We also thank the comments from František Gallovič, three other anonymous reviewers and associate editor, and Editor Rachel Abercrombie for checking the analysis process and improving the manuscript. Most of the computations are held at the French national computing center GENCI/CINES under the Grants A0050406700 (2018–2019), A0070406700 (2019–2020), and A0090406700 (2020–2021). A part of this work is funded by the “Institut National des Sciences de l’Univers” (IN-SU-CNRS), the “Réseau Sismologique & Géodésique Français” (RESIF), and the “Agence Nationale pour la Recherche” (ANR), under Grants ANR-11-EQPX-0040 and ANR-19-CE31-0003. S. Ruiz thanks the support of CONICYT/FONDECYT Project Number 1170430 and 1200779, and ProgramaRiesgoSísmico (Actividades de Interés Nacional [AIN]).
Barnhart, W. D., Murray, J. R. Briggs, R. W., Gomez, F., Miles, C. P. J., Svarc, J., et al. (2016). Coseismic slip and early afterslip of the 2015 Illapel, Chile, earthquake: Implications for frictional heterogeneity and coastal uplift. *Journal of Geophysical Research: Solid Earth*, **121**, 6172–6191. https://doi.org/10.1002/2016JB013124

Beck, S., Barrientos, S., Kausel, E., & Reyes, M. (1998). Source characteristics of historic earthquakes along the central Chile subduction Askew et alzone. *Journal of South American Earth Sciences*, **11**, 115–129. https://doi.org/10.1016/s0895-9811(98)00005-4

Bouchon, M., Karabulut, H., Aktar, M., Özalaybey, S., Schmittbuhl, J., Bouin, M.-P., & Marsan, D. (2021). The nucleation of the Izmit and Düzce earthquakes: Some mechanical logic on where and how ruptures began. *Geophysical Journal International*, **225**, 1510–1517. https://doi.org/10.1093/gji/ggab040

Carrasco, S., Ruiz, J. A., Contreras-Reyes, E., & Ortega-Culaciati, F. (2019). Shallow intraplate seismicity related to the Illapel 2015 Mw8.4 earthquake: Implications from the seismic source. *Tectonophysics*, **766**, 205–218. https://doi.org/10.1016/j.tecto.2019.06.011

Carvajal, M., Cisternas, M., & Catalán, P. A. (2017). Source of the 1730 Chilean earthquake from historical records: Implications for the future tsunami hazard on the coast of Metropolitan Chile. *Journal of Geophysical Research: Solid Earth*, **122**, 3648–3660. https://doi.org/10.1002/2017JB014063

Di Carli, S., François-Holden, C., Peyrat, S., & Madariaga, R. (2010). Dynamic inversion of the 2000 Tottori earthquake based on elliptical subfault approximations. *Journal of Geophysical Research*, **115**, B12238. https://doi.org/10.1029/2009JB006358

Di Giacomo, D., Engdahl, E. R., & Storchak, D. A. (2018). The ISC-GEM earthquake catalogue (1904–2014): Status after the extension project. *Earth System Science Data*, **10**, 1877–1899. https://doi.org/10.5194/essd-10-1877-2018

Dziewonski, A. M., Chou, T.-A., & Woodhouse, J. H. (1981). Determination of earthquake source parameters from waveform data for studies of global and regional seismicity. *Journal of Geophysical Research*, **86**, 2825–2852. https://doi.org/10.1029/JB086IB04p02825

Ekström, G., Nettles, M., & Dziewoński, A. M. (2012). The global CMT project 2004–2010: Centroid-moment tensors for 13,017 earthquakes. *Physics of the Earth and Planetary Interiors*, **200–201**, 1–9. https://doi.org/10.1016/j.pepi.2012.04.002

Ellsworth, W. L., & Beroza, G. C. (1995). Seismic evidence for an earthquake nucleation phase. *Science*, **268**, 851–855. https://doi.org/10.1126/science.268.5212.851

Ellsworth, W. L., & Bulut, F. (2018). Nucleation of the 1999 Izmit earthquake by a triggered cascade of foreshocks. *Nature Geoscience*, **11**, 531–535. https://doi.org/10.1038/s41561-018-0145-1

Fernández, J., Pastén, C., Ruiz, S., & Leyton, F. (2019). Damage assessment of the 2015 Mw8.3 Illapel earthquake in the north-central chile. *Natural Hazards*, **96**, 269–283. https://doi.org/10.1007/s11069-018-3541-3

Fuentes, M., Riquelme, S., Hayes, G., Medina, M., Melgar, D., Vargas, G., et al. (2017). A Study of the 2015 Mw8.3 Illapel earthquake and tsunami: Numerical and analytical approaches. In C. Braitenberg, & A. Rabinovich (Eds.), *The Chile-2015 (Illapel) earthquake and tsunami*. Birkhäuser, Cham: Pageoph. https://doi.org/10.1007/978-3-319-57822-4_18

Fukuyama, E., & Mikumo, T. (2007). Slip-weakening distance estimated at near-fault stations. *Geophysical Research Letters*, **34**, L09302. https://doi.org/10.1029/2006GL029203

Gallovič, F., Zahradník, J., Plicka, V., Sokos, E., Evangelidis, C., Foutoulakis, I., & Turhan, F. (2020). Complex rupture dynamics on an immature fault during the 2020 Mw6.8 Elazığ earthquake, Turkey. *Communications Earth and Environment*, **1**, 40. https://doi.org/10.1038/s43247-020-00038-x

Galvez, P., Ampuero, J.-P., Dalguer, L. A., Somala, S. N., & Nissen-Meyer, T. (2014). Dynamic earthquake rupture modelled with an un-structured 3-D spectral element method applied to the 2011 M9 Tohoku earthquake. *Geophysical Journal International*, **198**, 1222–1240. https://doi.org/10.1093/gji/ggu203

Goldberg, D. (1989). *Genetic algorithms in search, optimization, and machine learning* (372 pp.). Addison-Wesley Professional. ISBN978-0201157673.

Goldstein, P., & Snoke, A. (2005). *SAC availability for the IRIS community*. Incorporated Institutions for Seismology Data Management Center Electronic Newletter.

Grandin, R., Klein, E., Métois, M., & Vigny, C. (2016). Three-dimensional displacement field of the 2015 Mw8.3 Illapel earthquake (Chile) from across- and along-track Sentinel-1 TOPS interferometry. *Geophysical Research Letters*, **43**, 2552–2561. https://doi.org/10.1002/2016GL067954

Guatteri, M., & Spudich, P. (2000). What can strong-motion data tell us about slip-weakening fault-friction laws? *Bulletin of the Seismological Society of America*, **90**, 98–116. https://doi.org/10.1785/0119990053

Heidarzadeh, M., Murotani, S., Satake, K., Ishibe, T., & Gusman, A. R. (2016). Source model of the 16 September 2015 Illapel, Chile, Mw8.4 earthquake based on teleseismic and tsunami data. *Geophysical Research Letters*, **43**, 643–650. https://doi.org/10.1002/2015GL0672797

Henry, C., Das, S., & Woodhouse, J. H. (2000). The great March 25, 1998, Antarctic Plate earthquake: Moment tensor and rupture history. *Journal of Geophysical Research*, **105**, 16097–16118. https://doi.org/10.1029/2000JB900077

Hjörleifsdóttir, V., Kanamori, H., & Tromp, J. (2009). Modeling 3-D wave propagation and finite slip for the 1998 Balleny Islands earthquake. *Journal of Geophysical Research*, **114**, B03301. https://doi.org/10.1029/2008JB005975

Ida, Y. (1972). Cohesive force across the tip of a longitudinal-shear crack and Griffith's specific surface energy. *Journal of Geophysical Research*, **77**, 3796–3805. https://doi.org/10.1029/JB077i020p03796

Ide, S., & Aochi, H. (2005). Earthquakes as multiscale dynamic ruptures with heterogeneous fracture surface energy. *Journal of Geophysical Research*, **110**, B11303. https://doi.org/10.1029/2004JB003591

Ide, S., & Aochi, H. (2013). Historical seismicity and dynamic rupture process of the 2011 Tohoku-Oki earthquake. *Tectonophysics*, **600**, 1–13. https://doi.org/10.1016/j.tecto.2012.10.018

International Seismological Centre. (2018). ISC-GEM earthquake catalogue (v5.1). https://doi.org/10.31905/d808b825

Irikura, K., & Miyake, H. (2011). Recipe for predicting strong ground motion from crustal earthquake scenarios. *Pure and Applied Geophysics*, **168**, 85–104. https://doi.org/10.1007/s00024-010-0150-9

Kikuchi, M., & Kanamori, H. (1991). Inversion of complex body waves—III. *Bulletin of the Seismological Society of America*, **81**, 2335–2350.

Klein, E., Vigny, C., Fleitout, L., Grandin, R., Jolivet, R., Rivera, E., & Métois, M. (2017). A comprehensive analysis of the Illapel 2015 Mw8.3 earthquake from GPS and InSAR data. *Earth and Planetary Science Letters*, **469**, 123–134. https://doi.org/10.1016/j.epsl.2017.04.010

Kurahashi, S., & Irikura, K. (2011). Source model for generating strong ground motions during the 2011 off the pacific coast of Tohoku earthquake. *Earth, Planets and Space*, **63**, 571–576. https://doi.org/10.5047/eps.2011.06.044

Lange, D., Geersen, J., Barrientos, S., Moreno, M., Grevemeyer, I., Contreras-Reyes, E., & Kopp, H. (2016). Aftershock seismicity and tectonic setting of the 2015 September 16 Mw8.3 Illapel earthquake, Central Chile. *Geophysical Journal International*, **206**, 1424–1430. https://doi.org/10.1093/gji/ggw218

Lay, T., Li, L., & Cheung, K. F. (2016). Modeling tsunami observations to evaluate a proposed late tsunami earthquake stage for the 16 September 2015 Illapel, Chile, Mw8.3 earthquake. *Geophysical Research Letters*, **43**, 7902–7912. https://doi.org/10.1002/2016GL070002
Lemoine, A., Madariaga, R., & Campos, J. (2001). Evidence for earthquake interaction in central Chile: The July 1997-September 1998 sequence. *Geophysical Research Letters*, *28*, 2743-2746. https://doi.org/10.1029/2000GL012314

Leyton, F., Leopold, A., Hurtado, G., Pastén, C., Ruiz, S., Montalva, G., & Saéz, E. (2018a). Geophysical characterization of the Chilean seismological stations: First results. *Seismological Research Letters*, *89*, 519-525. https://doi.org/10.1785/0220170156

Leyton, F., Ruiz, S., Baez, J. C., Meneses, G., & Madariaga, R. (2018b). How fast can we reliably estimate the magnitude of subduction earthquakes? *Geophysical Research Letters*, *45*, 9633-9641. https://doi.org/10.1029/2018GL078991

Li, L., Lay, T., Cheung, K. F., & Ye, L. (2016). Joint modeling of teleseismic and tsunami wave observations to constrain the 16 September 2015 Illapel, Chile, Mw8.3 earthquake rupture process. *Geophysical Research Letters*, *43*, 4303-4312. https://doi.org/10.1002/2016GL068674

Madariaga, R., & Olsen, K. B. (2000). Criticality of rupture dynamics in 3-D. *Pure and Applied Geophysics*, *157*, 1981-2001. https://doi.org/10.1007/pl000001071

Mai, P. M., Spudich, P., & Boatwright, J. (2005). Hypocenter locations in finite-source rupture models. *Bulletin of the Seismological Society of America*, *95*, 965-980. https://doi.org/10.1785/0120040111

Malgrange, M., Deschamps, A., & Madariaga, R. (1981). Thrust and extensional faulting under the Chilean coast: 1965, 1971 Aconcagua earthquakes. *Geophysical Journal International*, *66*, 313-331. https://doi.org/10.1111/j.1365-246x.1981.tb05958.x

Melgar, D., Fan, W., Riquelme, S., Geng, J., Liang, C., Fuentes, M., et al. (2016). Slip segmentation and slow rupture to the trench during the 2015, Mw8.3 Illapel, Chile earthquake. *Geophysical Research Letters*, *43*, 961-966. https://doi.org/10.1002/2015GL067369

Meng, L., Bao, H., Huang, H., Zhang, A., Bloore, A., & Liu, Z. (2018). Double pincer movement: Encircling rupture splitting during the 2015 Mw8.3 Illapel earthquake. *Earth and Planetary Science Letters*, *495*, 164-173. https://doi.org/10.1016/j.epsl.2018.04.057

Métois, M., Vigny, C., & Socquet, A. (2016). Interseismic coupling, megathrust earthquakes and seismic swarms along the Chilean subduction zone (38°-18°S). *Pure and Applied Geophysics*, *173*, 1431-1449. https://doi.org/10.1007/s00024-016-1280-5

Ohnaka, M. (2003). A constitutive scaling law and a unified comprehension for frictional slip failure, shear fracture of intact rock, and earthquake rupture. *Journal of Geophysical Research*, *108*(B2), 2080. https://doi.org/10.1029/2000JB000123

Okuwaki, R., Yagi, Y., Aranguiz, R., González, J., & González, G. (2016). Rupture process during the 2015 Illapel, Chile earthquake: Zigzag-along-dip rupture episodes. In *Chile-2015 (Illapel) earthquake and tsunami*. (Vol. 173, pp. 1011-1020). Pageoph. https://doi.org/10.1007/s00024-016-1271-6

Pardo, M., Comte, D., Monfret, T., Boroschek, R., & Astroza, M. (2002). The October 15, 1997 Punitaqui earthquake (Mw=7.1): A destructive event within the subducting Nazca plate in central Chile. *Tectonophysics*, *345*, 199-210. https://doi.org/10.1016/s0040-1951(01)00213-x

Peyrat, S., Madariaga, R., Buforn, E., Campos, J., Asch, G., & Vilotte, J. P. (2010). Kinematic rupture process of the 2007 Tocopilla earthquake and its main aftershocks from teleseismic and strong-motion data. *Geophysical Journal International*, *182*, 1411-1430. https://doi.org/10.1111/j.1365-246X.2010.04685.x

Poli, P., Maksymowicz, A., & Ruiz, S. (2017). The Mw8.3 Illapel earthquake (Chile): Preseismic and postseismic activity associated with hydrated slab structures. *Geology*, *45*, 247-250. https://doi.org/10.1130/G38522.1

Ruiz, S., Aden-Antoniow, F., Baez, J. C., Otarola, C., Potin, B., del Campo, F., et al. (2017). Nucleation phase and dynamic inversion of the Mw6.9 Valparaiso 2017 earthquake in central Chile. *Geophysical Research Letters*, *44*, 10290-10297. https://doi.org/10.1002/2017GL075675

Ruiz, S., Klein, E., del Campo, F., Rivera, E., Poli, P., Metois, M., et al. (2016). The seismic sequence of the 16 September 2015 Mw8.3 Illapel, Chile, Earthquake. *Seismological Research Letters*, *87*, 789-799. https://doi.org/10.1785/0220150281

Ruiz, S., & Madariaga, R. (2011). Determination of the friction law parameters of the Mw6.7 Michilla earthquake in northern Chile by dynamic inversion. *Geophysical Research Letters*, *38*, L09317. https://doi.org/10.1029/2011GL047147

Ruiz, S., & Madariaga, R. (2013). Kinematic and dynamic inversion of the 2008 Northern Iwate earthquake. *Bulletin of the Seismological Society of America*, *103*, 694-708. https://doi.org/10.1785/010120056

Ruiz, S., & Madariaga, R. (2018). Historical and recent large megathrust earthquakes in Chile. *Tectonophysics*, *733*, 37-56. https://doi.org/10.1016/j.tecto.2018.01.015

Ruiz, S., Madariaga, R., Astroza, M., Saragoni, G. R., Lancieri, M., Vigny, C., & Campos, J. (2012). Short period rupture process of the 2010 Mw8.8 Maule earthquake in Chile. *Earthquake Spectra*, *28*(Supp.1), S1-S18. https://doi.org/10.1193/1.4000039

Saito, M. (1978). An automatic design algorithm for band selective recursive digital filters. *Butsuri Tanko (Journal of the Society of Exploration Geophysicists of Japan)*, *31*(4), 112-134 (in Japanese with English abstract)

Sánchez-Reyes, H. S., Tago, J., Métivier, L., Cruz-Atienza, V. M., & Virieux, J. (2018). An evolutive linear kinematic source inversion. *Journal of Geophysical Research: Solid Earth*, *123*, 4859-4885. https://doi.org/10.1029/2017/JB015388

Satake, K., & Heidarzadeh, M. (2017). A review of source models of the 2015 Illapel, Chile earthquake and insights from tsunami data. In C. Braitenberg, A. Rabinovich (Eds.). *The Chile-2015 (Illapel) Earthquake and Tsunami* (pp. 1-9). Birkhäuser, Cham: Pageoph. https://doi.org/10.1007/978-3-319-57822-4_1

Tilmann, F., Zhang, Y., Moreno, M., Saul, J., Eckelmann, F., Palo, M., et al. (2016). The 2015 Illapel earthquake, central Chile: A type case for a characteristic earthquake? *Geophysical Research Letters*, *43*, 574-583. https://doi.org/10.1002/2015GL066963

Tsai, V. C., Nettles, M., Ekstrom, G., & Dziewonski, A. M. (2005). Multiple CMT source analysis of the 2004 Sumatra earthquake. *Geophysical Research Letters*, *32*, L17304. https://doi.org/10.1029/2005GL023813

Twardzik, C., Das, S., & Madariaga, R. (2014). Inversion for the physical parameters that control the source dynamics of the 2004 Parkfield earthquake. *Journal of Geophysical Research: Solid Earth*, *119*, 7010-7027. https://doi.org/10.1002/2014JB011238

Uchide, T., & Ide, S. (2007). Development of multiscale slip inversion method and its application to the 2004 mid-Niigata Prefecture earthquake. *Journal of Geophysical Research*, *112*, B06313. https://doi.org/10.1029/2006JB004528

Udias, A., Madariaga, R., Buforn, E., Muñoz, D., & Ros, M. (2012). The large Chilean historical earthquakes of 1647, 1657, 1730, and 1751 from contemporary documents. *Bulletin of the Seismological Society of America*, *102*, 1639-1653. https://doi.org/10.1785/0120110289

Ulrich, T., & Aochi, H. (2015). Rapidness and robustness of finite-source inversion of the 2011 Mw9.0 Tohoku Earthquake by an elliptical-patches method using continuous GPS and acceleration data. *Pure and Applied Geophysics*, *172*, 3439-3453. https://doi.org/10.1007/s00024-014-0857-0

Ulrich, T., Gabriel, A. A., Ampuero, J. P., & Xu, W. (2019). Dynamic viability of the 2016 Mw7.8 Kaikōira earthquake cascade on weak crustal faults. *Nature Communications*, *10*, 1213. https://doi.org/10.1038/s41467-019-09125-w

Vallée, M., & Bouchon, M. (2004). Imaging coseismic rupture in far field by slip patches. *Geophysical Journal International*, *156*, 615-630. https://doi.org/10.1111/j.1365-246X.2004.02158.x

Voisin, C. (2002). Dynamic triggering of earthquakes: The nonlinear slip-dependent friction case. *Journal of Geophysical Research*, *107*, 2356. https://doi.org/10.1029/2001JB001121

Wang, R., Martin, F. L., & Roth, F. (2003). Computation of deformation induced by earthquakes in a multi-layered elastic crust-FORTRAN programs EDGRN/EDCMP. *Computers and Geosciences*, *29*, 195-207. https://doi.org/10.1016/S0098-3004(02)00111-5
Weng, H., & Yang, H. (2018). Constraining frictional properties on fault by dynamic rupture simulations and near-field observations. *Journal of Geophysical Research: Solid Earth*, **12**, 6658–6670. https://doi.org/10.1029/2017JB015414

Williamson, A., Newman, A., & Cummins, P. (2017). Reconstruction of coseismic slip from the 2015 Illapel earthquake using combined geodetic and tsunami waveform data. *Journal of Geophysical Research: Solid Earth*, **122**, 2119–2130. https://doi.org/10.1002/2016JB013883

Yao, S., & Yang, H. (2020). Rupture dynamics of the 2012 Nicoya Mw7.6 earthquake: Evidence for low strength on the megathrust. *Geophysical Research Letters*, **47**, e2020GL087508. https://doi.org/10.1029/2020GL087508

Ye, L., Lay, T., Kanamori, H., & Koper, K. D. (2016). Rapidly estimated seismic source parameters for the 16 September 2015 Illapel, Chile Mw8.3 Earthquake. *Pure and Applied Geophysics*, **173**, 321–332. https://doi.org/10.1007/s00024-015-1202-y

Yin, J., Yang, H., Yao, H., & Weng, H. (2016). Coseismic radiation and stress drop during the 2015 Mw8.3 Illapel, Chile megathrust earthquake. *Geophysical Research Letters*, **43**, 1520–1528. https://doi.org/10.1002/2015GL067381

Zahradník, J., Gallovic, F., Sokos, E., Serpetsidaki, A., & Tselentis, A. (2008). Quick fault-plane identification by a geometrical method: Application to the MW6.2 Leonidio earthquake, 6 January 2008, Greece. *Seismological Research Letters*, **79**, 653–662. https://doi.org/10.1785/gssrl.79.5.653

Zahradník, J., & Sokos, E. (2014). The Mw7.1 Van, Eastern Turkey, earthquake 2011: Two-point source modelling by iterative deconvolution and non-negative least squares. *Geophysical Journal International*, **196**, 522–538. https://doi.org/10.1093/gji/ggt386

Zhang, Y., Zhang, G., Hetland, E. A., Shan, X., Wen, S., & Zuo, R. (2016). Coseismic fault slip of the September 16, 2015 Mw8.3 Illapel, Chile earthquake estimated from InSAR data. *Pure and Applied Geophysics*, **173**, 1029–1038. https://doi.org/10.1007/s00024-016-1266-3