# Revisiting the 2015 Mw = 8.3 Illapel earthquake: unveiling complex fault slip properties using Bayesian inversion

E. Caballero, Z. Duputel, Cedric Twardzik, L. Rivera, E. Klein, J. Jiang, C. Liang, L. Zhu, R. Jolivet, E. Fielding, et al.

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

► To cite this version:

E. Caballero, Z. Duputel, Cedric Twardzik, L. Rivera, E. Klein, et al. Revisiting the 2015 Mw = 8.3 Illapel earthquake: unveiling complex fault slip properties using Bayesian inversion. Geophysical Journal International, 2023, 235, pp.2828-2845. 10.1093/gji/ggad380. insu-04582011

HAL Id: insu-04582011

https://insu.hal.science/insu-04582011

Submitted on 21 May 2024

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_2.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# Revisiting the 2015 Mw = 8.3 Illapel earthquake: unveiling complex fault slip properties using Bayesian inversion

E. Caballero, Z. Duputel, Cedric Twardzik, L. Rivera, E. Klein, J. Jiang, C. Liang, L. Zhu, R. Jolivet, E. Fielding, et al.

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

## To cite this version:

E. Caballero, Z. Duputel, Cedric Twardzik, L. Rivera, E. Klein, et al. Revisiting the 2015 Mw = 8.3 Illapel earthquake: unveiling complex fault slip properties using Bayesian inversion. Geophysical Journal International, 2023, 235 (3), pp.2828-2845. 10.1093/gji/ggad380. insu-04582011

HAL Id: insu-04582011

https://insu.hal.science/insu-04582011v1

Submitted on 21 May 2024

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_2.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
<figure><img src="image_3.png" /></figure>

# Revisiting the 2015 M<sub>w</sub> = 8.3 Illapel earthquake: unveiling complex fault slip properties using Bayesian inversion

E. Caballero<sup>1,2</sup> Z. Duputel<sup>1,4</sup> C. Twardzik,<sup>1,10</sup> L. Rivera<sup>1,5</sup> E. Klein,<sup>3</sup> J. Jiang<sup>1,5</sup> C. Liang,<sup>6</sup> L. Zhu<sup>1,7</sup> R. Jolivet,<sup>3,8</sup> E. Fielding<sup>9</sup> and M. Simons<sup>7</sup>

Accepted 2023 September 22. Received 2023 September 8; in original form 2023 March 31

## SUMMARY

The 2015 moment magnitude $M_w = 8.3$ Illapel earthquake is the largest mega-thrust earthquake that has been recorded along the Chilean subduction zone since the 2010 $M_w = 8.8$ Maule earthquake. Previous studies indicate a rupture propagation from the hypocentre to shallower parts of the fault, with a maximum slip varying from 10 to 16 m. The amount of shallow slip differs dramatically between rupture models with some results showing almost no slip at the trench and other models with significant slip at shallow depth. In this work, we revisit this event by combining a comprehensive data set including continuous and survey GNSS data corrected for post-seismic and aftershock signals, ascending and descending INSAR images of the Sentinel-1A satellite, tsunami data along with high-rate GPS, and doubly integrated strong-motion waveforms. We follow a Bayesian approach, in which the solution is an ensemble of models. The kinematic inversion is done using the cascading capability of the AlTar algorithm, allowing us to first get a static solution before integrating seismic data in a joint model. In addition, we explore a new approach to account for forward problem uncertainties using a second-order perturbation approach. Results show a rupture with two main slip patches, with significant slip at shallow depth. During the rupture propagation, we observe two regions that are encircled by the rupture, with no significant slip, westward of the hypocentre. These encircling effects have been previously suggested by back-projection results but have not been observed in finite-fault slip models. We propose that the encircled regions correspond to zones where the yield stress largely exceeds the initial stress or where fracture energy is too large to be ruptured during the Illapel earthquake. These asperities may potentially break in the future and probably already broke in the past.

**Key words:** Inverse theory; Probability distributions; Earthquake source observations.

## 1 INTRODUCTION

Chile is one of the most seismically active regions on Earth, where the Nazca Plate subducts under the South American Plate with a convergence rate of approximately 67 mm yr⁻<fnref n="1" /> (Angermann et al. 1999; Vigny et al. 2009). This large plate convergence rate is accommodated in parts by the occurrence of large megathrust earthquakes, such as the 1943 moment magnitude $M_w = 7.9–8.3$ Illapel event,

the 1960 $M_w = 9.5$ Valdivia earthquake, the 2010 $M_w = 8.8$ Maule earthquake, and the 2014 $M_w = 8.1$ Iquique earthquake (Lomnitz 2004; Ruiz & Madariaga 2018). The latest megathrust earthquake in Chile is the 2015 $M_w = 8.3$ Illapel earthquake, which occurred off the west coast of the Coquimbo region on 16 September 2015, at 22:54:31 UTC (Centro Sismológico Nacional, CSN; Li et al. 2016; Ruiz & Madariaga 2018). The 2015 Illapel earthquake initiated at a depth of 23 km and triggered a trans-pacific tsunami with

<footnote n="1"><sup>1</sup>Université de Strasbourg/EOST, CNRS, Institut Terre et Environment de Strasbourg, UMR 7063, 67084 Strasbourg, France.</footnote>

<footnote n="2"><sup>2</sup>Université de Grenoble-Alpes, CNRS, Institut des Sciences de la Terre, UMR 5275, 38000 Grenoble, France, Email: emmanuel.caballero-leyva@univ-grenoble-alpes.fr</footnote>

<footnote n="3"><sup>3</sup>Laboratoire de Géologie - CNRS UMR 8538, École Normale Supérieure - PSL University, 75230 Paris, France</footnote>

<footnote n="4"><sup>4</sup>Observatoire Volcanologique du Piton de La Fournaise, Université Paris Cité, Institut de Physique du Globe de Paris, CNRS, 97418 La Réunion, France</footnote>

<footnote n="5"><sup>5</sup>University of Oklahoma Norman Campus, School of Geosciences, Norman, OK 73019, United States</footnote>

<footnote n="6"><sup>6</sup>Institute of Remote Sensing and Geographic Information System, Peking University, Beijing 100871, China</footnote>

<footnote n="7"><sup>7</sup>Seismological Laboratory, California Institute of Technology, Pasadena, CA 91125, United States</footnote>

<footnote n="8"><sup>8</sup>Institut Universitaire de France, 1 rue Descartes, 75005 Paris, France</footnote>

<footnote n="9"><sup>9</sup>Jet Propulsion Laboratory, California Institute of Technology, Pasadena, CA 91109, United States</footnote>

<footnote n="10"><sup>10</sup>Université Côte d'Azur, CNRS, Observatoire de la Côte d'Azur, IRD, Geoazur, UMR 7329, 06560 Valbonne, France</footnote>
<figure><img src="image_3.png" /></figure>

# Revisiting the 2015 M<sub>w</sub> = 8.3 Illapel earthquake: unveiling complex fault slip properties using Bayesian inversion

E. Caballero<sup>1,2</sup> Z. Duputel<sup>1,4</sup> C. Twardzik,<sup>1,10</sup> L. Rivera<sup>1,5</sup> E. Klein,<sup>3</sup> J. Jiang<sup>1,5</sup> C. Liang,<sup>6</sup> L. Zhu<sup>1,7</sup> R. Jolivet,<sup>3,8</sup> E. Fielding<sup>9</sup> and M. Simons<sup>7</sup>

Accepted 2023 September 22. Received 2023 September 8; in original form 2023 March 31

## SUMMARY

The 2015 moment magnitude $M_w = 8.3$ Illapel earthquake is the largest mega-thrust earthquake that has been recorded along the Chilean subduction zone since the 2010 $M_w = 8.8$ Maule earthquake. Previous studies indicate a rupture propagation from the hypocentre to shallower parts of the fault, with a maximum slip varying from 10 to 16 m. The amount of shallow slip differs dramatically between rupture models with some results showing almost no slip at the trench and other models with significant slip at shallow depth. In this work, we revisit this event by combining a comprehensive data set including continuous and survey GNSS data corrected for post-seismic and aftershock signals, ascending and descending INSAR images of the Sentinel-1A satellite, tsunami data along with high-rate GPS, and doubly integrated strong-motion waveforms. We follow a Bayesian approach, in which the solution is an ensemble of models. The kinematic inversion is done using the cascading capability of the AlTar algorithm, allowing us to first get a static solution before integrating seismic data in a joint model. In addition, we explore a new approach to account for forward problem uncertainties using a second-order perturbation approach. Results show a rupture with two main slip patches, with significant slip at shallow depth. During the rupture propagation, we observe two regions that are encircled by the rupture, with no significant slip, westward of the hypocentre. These encircling effects have been previously suggested by back-projection results but have not been observed in finite-fault slip models. We propose that the encircled regions correspond to zones where the yield stress largely exceeds the initial stress or where fracture energy is too large to be ruptured during the Illapel earthquake. These asperities may potentially break in the future and probably already broke in the past.

**Key words:** Inverse theory; Probability distributions; Earthquake source observations.

## 1 INTRODUCTION

Chile is one of the most seismically active regions on Earth, where the Nazca Plate subducts under the South American Plate with a convergence rate of approximately 67 mm yr⁻<fnref n="1" /> (Angermann et al. 1999; Vigny et al. 2009). This large plate convergence rate is accommodated in parts by the occurrence of large megathrust earthquakes, such as the 1943 moment magnitude $M_w = 7.9–8.3$ Illapel event,

the 1960 $M_w = 9.5$ Valdivia earthquake, the 2010 $M_w = 8.8$ Maule earthquake, and the 2014 $M_w = 8.1$ Iquique earthquake (Lomnitz 2004; Ruiz & Madariaga 2018). The latest megathrust earthquake in Chile is the 2015 $M_w = 8.3$ Illapel earthquake, which occurred off the west coast of the Coquimbo region on 16 September 2015, at 22:54:31 UTC (Centro Sismológico Nacional, CSN; Li et al. 2016; Ruiz & Madariaga 2018). The 2015 Illapel earthquake initiated at a depth of 23 km and triggered a trans-pacific tsunami with

<footnote n="1"><sup>1</sup>Université de Strasbourg/EOST, CNRS, Institut Terre et Environment de Strasbourg, UMR 7063, 67084 Strasbourg, France.</footnote>

<footnote n="2"><sup>2</sup>Université de Grenoble-Alpes, CNRS, Institut des Sciences de la Terre, UMR 5275, 38000 Grenoble, France, Email: emmanuel.caballero-leyva@univ-grenoble-alpes.fr</footnote>

<footnote n="3"><sup>3</sup>Laboratoire de Géologie - CNRS UMR 8538, École Normale Supérieure - PSL University, 75230 Paris, France</footnote>

<footnote n="4"><sup>4</sup>Observatoire Volcanologique du Piton de La Fournaise, Université Paris Cité, Institut de Physique du Globe de Paris, CNRS, 97418 La Réunion, France</footnote>

<footnote n="5"><sup>5</sup>University of Oklahoma Norman Campus, School of Geosciences, Norman, OK 73019, United States</footnote>

<footnote n="6"><sup>6</sup>Institute of Remote Sensing and Geographic Information System, Peking University, Beijing 100871, China</footnote>

<footnote n="7"><sup>7</sup>Seismological Laboratory, California Institute of Technology, Pasadena, CA 91125, United States</footnote>

<footnote n="8"><sup>8</sup>Institut Universitaire de France, 1 rue Descartes, 75005 Paris, France</footnote>

<footnote n="9"><sup>9</sup>Jet Propulsion Laboratory, California Institute of Technology, Pasadena, CA 91109, United States</footnote>

<footnote n="10"><sup>10</sup>Université Côte d'Azur, CNRS, Observatoire de la Côte d'Azur, IRD, Geoazur, UMR 7329, 06560 Valbonne, France</footnote>
waves reaching more than 4 m high in Chile (An & Meng 2017; Fernández et al. 2019). The thrust focal mechanism is consistent with the rupture of the megathrust interface (Ekström et al. 2012). Most source inversions suggested that the rupture lasted around 100 s (Heidarzadeh et al. 2016; Melgar et al. 2016; Tilmann et al. 2016) but some studies report much larger rupture durations (e.g. Lee et al. 2016). The previous earthquake to rupture this section of the megathrust occurred in 1943, with a smaller magnitude between $M_w = 7.9–8.3$, and a duration of approximately 30 s (Beck et al. 1998; Lomnitz 2004; Ruiz & Madariaga 2018). The hypocentral depth of the 1943 event is unfortunately not well resolved and is estimated between 10 and 30 km.

Different groups have published kinematic slip rupture models for the 2015 $M_w = 8.3$ Illapel earthquake. As discussed by Satake & Heidarzadeh (2017), even though all of these models share general features, some properties of the rupture are still under debate (Heidarzadeh et al. 2016; Li et al. 2016; Ruiz et al. 2016; Tilmann et al. 2016; Williamson et al. 2017; An & Meng 2017). For example, An & Meng (2017) suggest the absence of shallow slip, while other studies indicate that shallow slip is necessary to explain tsunami records (Lay et al. 2016; Li et al. 2016; Tilmann et al. 2016). In fact, Tilmann et al. (2016) suggested that the 1943 and 2015 events differ in their shallow slip.

The degree of rupture complexity also varies among previously published results. In contrast with the relatively simple rupture processes suggested by the aforementioned results, other studies suggest a more complex rupture scenario with at least two main slip asperities (Melgar et al. 2016; Lee et al. 2016). Although the relatively compact model of Melgar et al. (2016) is consistent with tsunami observations, Lay et al. (2016) show that the model of Lee et al. (2016) involving a broad area of shallow slip rupturing multiple times cannot reproduce tsunami data. Several back-projections studies confirm the complexity of the 2015 Illapel rupture (Melgar et al. 2016; Okuwaki et al. 2016; Yin et al. 2016). A common result among back-projection studies is that the Illapel earthquake presents a northwestward migration. For example, An et al. (2017) shows a complex frequency dependent rupture propagation with several branches. The back-projected low-frequency (LF) sources migrate mainly updip to the west, while the high-frequency (HF) sources initially move downdip towards the northeast before veering updip towards the northwest. On the other hand, Melgar et al. (2018) suggest a rupture that splits into two different branches separated along dip. The analysis of these multiple rupture branches suggests an encircling rupture that seems to be aligned with regions experiencing a high slip rate and large shallow slip. Unfortunately, such a complex pattern has not been confirmed by kinematic slip inversion models yet. Potentially, such encircling rupture effect is only constrained by the HF wavefield, hence not resolvable with slip inversions. In addition, such encircling pattern likely involves abrupt changes in rupture velocities, while most slip inversions consider fixed rupture velocities and smoothing constraints.

In this work, we revisit the 2015 $M_w = 8.3$ Illapel earthquake by combining a comprehensive data set including permanent and survey GPS stations corrected for post-seismic and aftershock signals, ascending and descending Sentinel-1A InSAR images along with high-rate GPS and doubly integrated strong-motion waveforms. We follow a Bayesian approach using the AlTar code, which allows us to obtain the posterior probability distribution of slip models rather than a single optimum solution. We also use a non-linear parametrization enabling significant variation of rupture velocity

during the rupture process. We also analyse the impact that prediction error covariance matrices have on coseismic slip inversions results.

## 2 DATA

We investigate the complex rupture of the 2015 $M_w = 8.3$ Illapel earthquake using multiple data sets that are shown in Fig. 1. This database includes GPS offsets, Interferometric Synthetic Aperture Radar (InSAR) images, tsunami data along with high-rate GPS and strong motion waveforms.

InSAR images are obtained from the Sentinel-1A satellite with ascending and descending orbits (see text S1). We use 14 tsunami stations: 6 DART buoys and 6 coast gauges focusing mainly on first arrivals and open sea sites to minimize coastal effects (see text S2). We use daily and survey GPS data provided by Klein et al. (2017). Both data sets are affected by co-seismic offsets induced by $M_w = 7.1$ and $M_w = 6.8$ aftershocks occurring, respectively, 23 min and 5 hr after the main shock. Survey GPS data also includes several weeks of post-seismic displacement. Details of GPS data processing can be found in Klein et al. (2017). To correct both daily and campaign GPS data from aftershocks and post-seismic deformation, we use high-rate post-seismic time-series from Twardzik et al. (2021). These measurements are spatially interpolated using cubic splines and removed from coseismic GPS offsets. We estimate uncertainty associated with the corrected data by conducting the aforementioned correction stochastically (using Gaussian realizations given uncertainties on daily, survey and post-seismic GPS data sets). A comparison between corrected and uncorrected GPS data is shown in Fig. A1. We note that the nominal standard deviations of the GPS data are unrealistically small (i.e. on the order of 5–10 mm), leading to overfitting of the GPS coseismic displacements in the inversion procedure. To mitigate this issue, we scale the resulting standard deviations to ensure a unit reduced $\chi^2$, a statistical indicator that helps to correct for over or under estimation of uncertainties (Supporting Information text S3). As a result, we increase the standard deviation of the GPS static displacements by a factor of 10 for the East component and 5 for the north and vertical components. Although this approach is empirical, it allows us to avoid any overfitting of the GPS observations while keeping a relative weighting between stations based on the variability of the corrected observations.

For the kinematic data set (i.e. seismic waveforms), we use records from High Rate GPS (HRGPS) stations and strong motion data located within 5° from the main shock hypocentre. These stations are part of the Chilean Seismological Service (CSN) of the Universidad de Chile (Universidad de Chile 2012). In total, we have 96 strong motion waveforms that we double integrate into displacement time series and 12 HRGPS components. The integration of acceleration data is a delicate operation that can easily result in large drifts in velocity and displacement waveforms. Therefore, to obtain displacement records, after removing any linear trend in accelerograms, we remove an additional velocity drift at the end of the waveforms. This additional coda correction is done by using a quadratic function to fit displacement waveforms from the time when 90 per cent of the acceleration energy is reached. Visual inspection of the corrected displacement records is then done to ensure the good quality of the data. To further check the corrected records, we compare the obtained strong motion displacements with HRGPS displacements (Figs 2 and A2). In total, we were able to recover 43 displacement components from strong motion with high-quality displacement waveforms.
waves reaching more than 4 m high in Chile (An & Meng 2017; Fernández et al. 2019). The thrust focal mechanism is consistent with the rupture of the megathrust interface (Ekström et al. 2012). Most source inversions suggested that the rupture lasted around 100 s (Heidarzadeh et al. 2016; Melgar et al. 2016; Tilmann et al. 2016) but some studies report much larger rupture durations (e.g. Lee et al. 2016). The previous earthquake to rupture this section of the megathrust occurred in 1943, with a smaller magnitude between $M_w = 7.9–8.3$, and a duration of approximately 30 s (Beck et al. 1998; Lomnitz 2004; Ruiz & Madariaga 2018). The hypocentral depth of the 1943 event is unfortunately not well resolved and is estimated between 10 and 30 km.

Different groups have published kinematic slip rupture models for the 2015 $M_w = 8.3$ Illapel earthquake. As discussed by Satake & Heidarzadeh (2017), even though all of these models share general features, some properties of the rupture are still under debate (Heidarzadeh et al. 2016; Li et al. 2016; Ruiz et al. 2016; Tilmann et al. 2016; Williamson et al. 2017; An & Meng 2017). For example, An & Meng (2017) suggest the absence of shallow slip, while other studies indicate that shallow slip is necessary to explain tsunami records (Lay et al. 2016; Li et al. 2016; Tilmann et al. 2016). In fact, Tilmann et al. (2016) suggested that the 1943 and 2015 events differ in their shallow slip.

The degree of rupture complexity also varies among previously published results. In contrast with the relatively simple rupture processes suggested by the aforementioned results, other studies suggest a more complex rupture scenario with at least two main slip asperities (Melgar et al. 2016; Lee et al. 2016). Although the relatively compact model of Melgar et al. (2016) is consistent with tsunami observations, Lay et al. (2016) show that the model of Lee et al. (2016) involving a broad area of shallow slip rupturing multiple times cannot reproduce tsunami data. Several back-projections studies confirm the complexity of the 2015 Illapel rupture (Melgar et al. 2016; Okuwaki et al. 2016; Yin et al. 2016). A common result among back-projection studies is that the Illapel earthquake presents a northwestward migration. For example, An et al. (2017) shows a complex frequency dependent rupture propagation with several branches. The back-projected low-frequency (LF) sources migrate mainly updip to the west, while the high-frequency (HF) sources initially move downdip towards the northeast before veering updip towards the northwest. On the other hand, Melgar et al. (2018) suggest a rupture that splits into two different branches separated along dip. The analysis of these multiple rupture branches suggests an encircling rupture that seems to be aligned with regions experiencing a high slip rate and large shallow slip. Unfortunately, such a complex pattern has not been confirmed by kinematic slip inversion models yet. Potentially, such encircling rupture effect is only constrained by the HF wavefield, hence not resolvable with slip inversions. In addition, such encircling pattern likely involves abrupt changes in rupture velocities, while most slip inversions consider fixed rupture velocities and smoothing constraints.

In this work, we revisit the 2015 $M_w = 8.3$ Illapel earthquake by combining a comprehensive data set including permanent and survey GPS stations corrected for post-seismic and aftershock signals, ascending and descending Sentinel-1A InSAR images along with high-rate GPS and doubly integrated strong-motion waveforms. We follow a Bayesian approach using the AlTar code, which allows us to obtain the posterior probability distribution of slip models rather than a single optimum solution. We also use a non-linear parametrization enabling significant variation of rupture velocity

during the rupture process. We also analyse the impact that prediction error covariance matrices have on coseismic slip inversions results.

## 2 DATA

We investigate the complex rupture of the 2015 $M_w = 8.3$ Illapel earthquake using multiple data sets that are shown in Fig. 1. This database includes GPS offsets, Interferometric Synthetic Aperture Radar (InSAR) images, tsunami data along with high-rate GPS and strong motion waveforms.

InSAR images are obtained from the Sentinel-1A satellite with ascending and descending orbits (see text S1). We use 14 tsunami stations: 6 DART buoys and 6 coast gauges focusing mainly on first arrivals and open sea sites to minimize coastal effects (see text S2). We use daily and survey GPS data provided by Klein et al. (2017). Both data sets are affected by co-seismic offsets induced by $M_w = 7.1$ and $M_w = 6.8$ aftershocks occurring, respectively, 23 min and 5 hr after the main shock. Survey GPS data also includes several weeks of post-seismic displacement. Details of GPS data processing can be found in Klein et al. (2017). To correct both daily and campaign GPS data from aftershocks and post-seismic deformation, we use high-rate post-seismic time-series from Twardzik et al. (2021). These measurements are spatially interpolated using cubic splines and removed from coseismic GPS offsets. We estimate uncertainty associated with the corrected data by conducting the aforementioned correction stochastically (using Gaussian realizations given uncertainties on daily, survey and post-seismic GPS data sets). A comparison between corrected and uncorrected GPS data is shown in Fig. A1. We note that the nominal standard deviations of the GPS data are unrealistically small (i.e. on the order of 5–10 mm), leading to overfitting of the GPS coseismic displacements in the inversion procedure. To mitigate this issue, we scale the resulting standard deviations to ensure a unit reduced $\chi^2$, a statistical indicator that helps to correct for over or under estimation of uncertainties (Supporting Information text S3). As a result, we increase the standard deviation of the GPS static displacements by a factor of 10 for the East component and 5 for the north and vertical components. Although this approach is empirical, it allows us to avoid any overfitting of the GPS observations while keeping a relative weighting between stations based on the variability of the corrected observations.

For the kinematic data set (i.e. seismic waveforms), we use records from High Rate GPS (HRGPS) stations and strong motion data located within 5° from the main shock hypocentre. These stations are part of the Chilean Seismological Service (CSN) of the Universidad de Chile (Universidad de Chile 2012). In total, we have 96 strong motion waveforms that we double integrate into displacement time series and 12 HRGPS components. The integration of acceleration data is a delicate operation that can easily result in large drifts in velocity and displacement waveforms. Therefore, to obtain displacement records, after removing any linear trend in accelerograms, we remove an additional velocity drift at the end of the waveforms. This additional coda correction is done by using a quadratic function to fit displacement waveforms from the time when 90 per cent of the acceleration energy is reached. Visual inspection of the corrected displacement records is then done to ensure the good quality of the data. To further check the corrected records, we compare the obtained strong motion displacements with HRGPS displacements (Figs 2 and A2). In total, we were able to recover 43 displacement components from strong motion with high-quality displacement waveforms.
<figure><img src="image_4.png" /><figcaption>**Figure 1.** General overview of the studied region with data sets used in this study (a). Green star represents the hypocentre obtained by the Chilean Seismological Center (CSN). White rectangles represent the fault geometry used in this study. Focal mechanisms correspond to aftershocks Global CMT solutions. Ascending (b) and descending (c) Sentinel-1A InSAR images. Small black arrows represent the LOS and orbit direction, respectively.</figcaption></figure>

To calculate synthetic static displacements, we use the Classic Slip Inversion (CSI) package (https://github.com/jolivet/csi), using the approach of Zhu & Rivera (2002) for a layered earth model. We calculate Green's functions using the 1-D velocity model built by Duputel et al. (2015) (see Fig. 3). For the kinematic Green's functions, we use the wavenumber integration code of the CPS seismology package (http://www.eas.slu.edu/eqc/eqccps.html) from Herrmann (2013). We filter both the kinematic Green's function and data in the 0.01–0.06667 Hz passband.

## 3 METHODOLOGY

To perform the inversion, we follow a Bayesian approach in which we obtain an ensemble of models and not a unique solution. The inversion is done using the cascading capability of the AlTar code (https://altar.readthedocs.io), allowing us to first get a static solution, and then to integrate waveform data in a joint model. This code is based on the Cascading Adaptative Metropolis In Parallel (CATMIP) algorithm proposed by Minson et al. (2013) that we will describe below. The AlTar package has been successfully used for different problems. Jolivet et al. (2015, 2020, 2023) estimated the interseismic coupling of the San Andreas fault, the Northern Chile

subduction interface and the North Anatolian fault. Studies of individual earthquakes have been carried out by Duputel et al. (2015), Bltery et al. (2016) and Gombert et al. (2018a), among others.

Starting from Bayes theorem, we write the a posteriori probability density function (PDF) of the parameters **m**, given the observations **d**obs:

$$p(\mathbf{m}|\mathbf{d}_{\text{obs}}) = \kappa p(\mathbf{m}) p(\mathbf{d}_{\text{obs}}|\mathbf{m}), \quad (1)$$

where $p(\mathbf{m})$ is the a priori probability density function of parameters, $p(\mathbf{d}_{\text{obs}}|\mathbf{m})$ is the data likelihood function and $\kappa$ a normalization factor. We define the likelihood function as:

$$p(\mathbf{d}_{\text{obs}}|\mathbf{m}) = \exp \left( -\frac{1}{2} (\mathbf{d}_{\text{obs}} - \mathbf{g}(\mathbf{m}))^{\mathrm{T}} \mathbf{C}_{\chi}^{-1} (\mathbf{d}_{\text{obs}} - \mathbf{g}(\mathbf{m})) \right). \quad (2)$$

$C_\chi$ is the misfit covariance matrix that is the sum of $C_d$ and $C_p$, which correspond to covariance matrices describing observational and forward modelling uncertainties, respectively. We sample the a posteriori PDF using a series of transitional intermediate PDF. The transitional PDFs are controlled by the tempering parameter $\beta$, which modulates the information content at each transitional step such as:

$$f(\mathbf{m}|\mathbf{d}_{\text{obs}}, \beta_k) = \kappa p(\mathbf{m}) p(\mathbf{d}_{\text{obs}}|\mathbf{m})^{\beta_k}, \quad (3)$$
<figure><img src="image_4.png" /><figcaption>**Figure 1.** General overview of the studied region with data sets used in this study (a). Green star represents the hypocentre obtained by the Chilean Seismological Center (CSN). White rectangles represent the fault geometry used in this study. Focal mechanisms correspond to aftershocks Global CMT solutions. Ascending (b) and descending (c) Sentinel-1A InSAR images. Small black arrows represent the LOS and orbit direction, respectively.</figcaption></figure>

To calculate synthetic static displacements, we use the Classic Slip Inversion (CSI) package (https://github.com/jolivet/csi), using the approach of Zhu & Rivera (2002) for a layered earth model. We calculate Green's functions using the 1-D velocity model built by Duputel et al. (2015) (see Fig. 3). For the kinematic Green's functions, we use the wavenumber integration code of the CPS seismology package (http://www.eas.slu.edu/eqc/eqccps.html) from Herrmann (2013). We filter both the kinematic Green's function and data in the 0.01–0.06667 Hz passband.

## 3 METHODOLOGY

To perform the inversion, we follow a Bayesian approach in which we obtain an ensemble of models and not a unique solution. The inversion is done using the cascading capability of the AlTar code (https://altar.readthedocs.io), allowing us to first get a static solution, and then to integrate waveform data in a joint model. This code is based on the Cascading Adaptative Metropolis In Parallel (CATMIP) algorithm proposed by Minson et al. (2013) that we will describe below. The AlTar package has been successfully used for different problems. Jolivet et al. (2015, 2020, 2023) estimated the interseismic coupling of the San Andreas fault, the Northern Chile

subduction interface and the North Anatolian fault. Studies of individual earthquakes have been carried out by Duputel et al. (2015), Bletery et al. (2016) and Gombert et al. (2018a), among others.

Starting from Bayes theorem, we write the a posteriori probability density function (PDF) of the parameters **m**, given the observations **d**<sub>obs</sub>:

$$p(\mathbf{m}|\mathbf{d}_{\text{obs}}) = \kappa p(\mathbf{m}) p(\mathbf{d}_{\text{obs}}|\mathbf{m}), \quad (1)$$

where $p(\mathbf{m})$ is the a priori probability density function of parameters, $p(\mathbf{d}_{\text{obs}}|\mathbf{m})$ is the data likelihood function and $\kappa$ a normalization factor. We define the likelihood function as:

$$p(\mathbf{d}_{\text{obs}}|\mathbf{m}) = \exp \left( -\frac{1}{2} (\mathbf{d}_{\text{obs}} - \mathbf{g}(\mathbf{m}))^{\mathrm{T}} \mathbf{C}_{\chi}^{-1} (\mathbf{d}_{\text{obs}} - \mathbf{g}(\mathbf{m})) \right). \quad (2)$$

$C_{\chi}$ is the misfit covariance matrix that is the sum of $C_d$ and $C_p$, which correspond to covariance matrices describing observational and forward modelling uncertainties, respectively. We sample the a posteriori PDF using a series of transitional intermediate PDF. The transitional PDFs are controlled by the tempering parameter $\beta$, which modulates the information content at each transitional step such as:

$$f(\mathbf{m}|\mathbf{d}_{\text{obs}}, \beta_k) = \kappa p(\mathbf{m}) p(\mathbf{d}_{\text{obs}}|\mathbf{m})^{\beta_k}, \quad (3)$$
<figure><img src="image_5.png" /><figcaption>**Figure 2.** Comparison between displacements corrected from strong motion records and HRGPS displacements. Red and black waveforms represent HRGPS and strong motion, respectively. On the maps, the blue star represents the CSN hypocentre while circles indicate station location (orange for the strong motion station considered, yellow for the other strong motion stations, and purple for HRGPS stations). φ, and Δ represent the azimuth and distance from the epicentre. The angle α is the component azimuth (0° north, 90° east). Time-shifts between waveforms are due to slight differences in station location (i.e. between HRGPS and strong motion records). Other examples of comparison are shown in Fig. A2.</figcaption></figure>

where (k = 1, ..., M) and β varies from zero to one, that is 0 = β₀
< β₁, ..., βₘ = 1.

These transitional steps will converge to the final solution by
smoothly informing the system (i.e. by increasing β). In addition,
we apply a cascading approach to improve the convergence of the
sampling by first solving for the static problem before sampling the
full joint kinematic slip inversion. More details about the algorithm
can be found in Minson *et al.* (2013). As mentioned before, the Cₓ
matrix incorporates different uncertainty assessments. The observa-
tional uncertainty is commonly related to errors in measurements.
The details of observational uncertainty estimates can be found in
text S4.

Prediction uncertainties are associated with imperfect forward
modelling that can be caused by different factors, such as imperfect
earth models or fault geometries (Wald & Graves 2001; Beres-
nev 2003; Ide 2015; Williams & Wallace 2015). Several studies
have highlighted the importance of considering forward modelling
uncertainties in slip inversions (Yagi & Fukahata 2011; Duputel
et al. 2012, 2014; Hallo & Gallovič 2016; Ragon et al. 2018). For
example, Duputel et al. (2014) study the uncertainties linked to in-
accuracies in the Earth structure model. On the other side, Ragon
et al. (2018) analyse uncertainties associated with inaccuracies in

fault geometries. Also, Razafindrakoto & Mai (2014) assess the in-
fluence of the employed source time function and elastic structure
on earthquake slip imaging.

In this study, we focus on accounting uncertainties due to Earth
structure modelling. Specifically, we evaluate the impact of inaccu-
racies in the 1-D velocity model employed to compute static and
kinematic predictions. Uncertainties in the elastic parameters Ψ is
assumed to follow a log-normal distribution:

$$
p(\log \Psi) = \frac{1}{\sqrt{(2\pi)^N |C_{\Psi}|}} \exp \left( -\frac{1}{2} (\log \Psi - \log \tilde{\Psi})^T C_{\Psi} (\log \Psi - \log \tilde{\Psi}) \right), \quad (4)
$$

where C_Ψ is the covariance characterizing uncertainty around log Ψ
(the logarithm of the elastic parameters used to compute the pre-
dictions shown in Fig. 3). This choice of a log-normal distribution
is motivated by the fact that (1) the elastic parameters are strictly
positive and (2) Ψ values are derived from tomography techniques
based on relative model perturbations (δlog Ψ; e.g. Tromp et al.
2005). The earth model uncertainty considered in this study is
shown in Fig. 3. This level of variability is measured by com-
paring different models from the region (following Duputel et al.
2015).
<figure><img src="image_5.png" /><figcaption>**Figure 2.** Comparison between displacements corrected from strong motion records and HRGPS displacements. Red and black waveforms represent HRGPS and strong motion, respectively. On the maps, the blue star represents the CSN hypocentre while circles indicate station location (orange for the strong motion station considered, yellow for the other strong motion stations, and purple for HRGPS stations). $\phi$, and $\Delta$ represent the azimuth and distance from the epicentre. The angle $\alpha$ is the component azimuth (0° north, 90° east). Time-shifts between waveforms are due to slight differences in station location (i.e. between HRGPS and strong motion records). Other examples of comparison are shown in Fig. A2.</figcaption></figure>

where ($k = 1, \dots, M$) and $\beta$ varies from zero to one, that is $0 = \beta_0 < \beta_1, \dots, \beta_M = 1$.

These transitional steps will converge to the final solution by smoothly informing the system (i.e. by increasing $\beta$). In addition, we apply a cascading approach to improve the convergence of the sampler by first solving for the static problem before sampling the full joint kinematic slip inversion. More details about the algorithm can be found in Minson et al. (2013). As mentioned before, the $C_\chi$ matrix incorporates different uncertainty assessments. The observational uncertainty is commonly related to errors in measurements. The details of observational uncertainty estimates can be found in text S4.

Prediction uncertainties are associated with imperfect forward modelling that can be caused by different factors, such as imperfect earth models or fault geometries (Wald & Graves 2001; Beresnev 2003; Ide 2015; Williams & Wallace 2015). Several studies have highlighted the importance of considering forward modelling uncertainties in slip inversions (Yagi & Fukahata 2011; Duputel et al. 2012, 2014; Hallo & Gallovič 2016; Ragon et al. 2018). For example, Duputel et al. (2014) study the uncertainties linked to inaccuracies in the Earth structure model. On the other side, Ragon et al. (2018) analyse uncertainties associated with inaccuracies in

fault geometries. Also, Razafindrakoto & Mai (2014) assess the influence of the employed source time function and elastic structure on earthquake slip imaging.

In this study, we focus on accounting uncertainties due to Earth structure modelling. Specifically, we evaluate the impact of inaccuracies in the 1-D velocity model employed to compute static and kinematic predictions. Uncertainties in the elastic parameters $\Psi$ is assumed to follow a log-normal distribution:

$$p(\log \Psi) = \frac{1}{\sqrt{(2\pi)^N |C_{\Psi}|}} \exp \left( -\frac{1}{2} (\log \Psi - \log \tilde{\Psi})^T C_{\Psi} (\log \Psi - \log \tilde{\Psi}) \right), \quad (4)$$

where $C_\Psi$ is the covariance characterizing uncertainty around $\log \tilde{\Psi}$ (the logarithm of the elastic parameters used to compute the predictions shown in Fig. 3). This choice of a log-normal distribution is motivated by the fact that (1) the elastic parameters are strictly positive and (2) $\Psi$ values are derived from tomography techniques based on relative model perturbations ($\delta \log \Psi$; e.g. Tromp et al. 2005). The earth model uncertainty considered in this study is shown in Fig. 3. This level of variability is measured by comparing different models from the region (following Duputel et al. 2015).
<figure><img src="image_6.png" /><figcaption>**Figure 3.** Model variability of the *P* wave, *S* wave and density as a function of depth in the Illapel region. The black line represents the velocity layered model used for Green's function (GF) calculation. Grey histograms are the probability density function for each parameter as a function of depth.</figcaption></figure>

<table><caption>**Table 1.** Approaches to calculate $C_p$ (for 36 parameters).</caption><thead><tr><th>Approach</th><th>Number of forward model evaluations</th></tr></thead><tbody><tr><td>Without C<sub>p</sub></td><td>0</td></tr><tr><td>Empirical</td><td>195 (in this study)</td></tr><tr><td>1st order forward derivatives</td><td>37</td></tr><tr><td>1st order centred derivatives</td><td>72</td></tr><tr><td>2nd order without cross-terms</td><td>73</td></tr><tr><td>2nd order</td><td>1333</td></tr></tbody></table>

We follow three different schemes to map earth model uncertainty into prediction uncertainty. The first straightforward approach is to empirically calculate the prediction uncertainty covariance matrix $C_p$ using predictions computed for a large number of random earth models $\Psi^i$, ($i = 1, \dots, n$) drawn from $p(\log \Psi)$:

$$C_p = \frac{1}{n-1} \sum_{i=1}^{n} (g(\Psi^i, \mathbf{m}) - g(\tilde{\Psi}, \mathbf{m})) (g(\Psi^i, \mathbf{m}) - g(\tilde{\Psi}, \mathbf{m}))^T, (5)$$

where $g(\Psi^i, \mathbf{m})$ is the prediction for the earth model $\Psi^i$ and the source model $\mathbf{m}$. In our case, we use a preliminary source model $\mathbf{m}$ derived from a first preliminary slip inversion. $g(\tilde{\Psi}, \mathbf{m})$ is the prediction response for the average earth model $\tilde{\Psi}$. This empirical approach is computationally expensive because it needs the calculation of predictions for each randomly generated earth model. To evaluate the number of models $n$ necessary to calculate an accurate empirical $C_p$ matrix, we compare empirical $C_p$ matrices calculated for an increasing number of random earth models. We observe that the empirical $C_p$ matrix is converging using 195 random Earth samples (Fig. A3), corresponding to relatively smooth histograms in Fig. 3.

To test a computationally less expensive approach, we also follow the first-order approximation approach proposed by Duputel *et al.* (2014). Assuming that we can approximate our forward model $g(\Psi, \mathbf{m})$ by linearized perturbations, for an *a priori* earth model we have

then:

$$g(\Psi, \mathbf{m}) \approx g(\tilde{\Psi}, \mathbf{m}) + K_{\Psi}(\tilde{\Psi}, \mathbf{m}) \cdot (\Psi - \tilde{\Psi}), \quad (6)$$

where $K$ is the sensitivity kernels of the prediction with respect to elastic parameters used to compute forward predictions:

$$(K_{\Psi})_{ij}(\tilde{\Psi}, \mathbf{m}) = \frac{\partial g_i}{\partial \Psi_j}(\tilde{\Psi}, \mathbf{m}), \quad (7)$$

where $\Psi_j$ corresponds to the jth elastic parameter in the earth model $\Psi$. We use then $K$ to estimate $C_p$ as:

$$C_p = K_{\Psi} \cdot C_{\Psi} \cdot K_{\Psi}^{T}, \quad (8)$$

where $C_{\Psi}$ is the same log-normal covariance that we use for perturbing the random models of the empirical $C_p$ in eq. (4). Although this approach looks appropriate for static data, it could be problematic for kinematic data as the link between earth model perturbations and waveform predictions is probably not linear. Indeed, changes in the velocity model induce both time-shifts and amplitude variations in the predicted waveforms.

Therefore, we also explore the possibility of using a 2nd order perturbation approach of the forward model as:

$$g(\Psi, \mathbf{m}) \approx g(\tilde{\Psi}, \mathbf{m}) + K_{\Psi}(\tilde{\Psi}, \mathbf{m}) \cdot (\Psi - \tilde{\Psi}) + \frac{1}{2} (\Psi - \tilde{\Psi}) \cdot H_{\Psi}(\tilde{\Psi}, \mathbf{m}) \cdot (\Psi - \tilde{\Psi}), \quad (9)$$

where $H_{\Psi}$ includes the second order derivative with respect to the elastic parameters:

$$(H_{\Psi})_{ijk}(\tilde{\Psi}, \mathbf{m}) = \frac{\partial^2 g_i}{\partial \Psi_k \partial \Psi_j}(\tilde{\Psi}, \mathbf{m}). \quad (10)$$

From eq. (9), we can then calculate the $C_p$ matrix using eq. (5) by rapidly generating a large number of forward model predictions.

The derivatives in eq. (9) are computed numerically using finite differences. We summarize the difference in computational cost between approaches in table 1. The computational cost of each approach in terms of forward model evaluation is summarized in Table 1. In this study, the empirical approach necessitated about 200 forward model evaluations, which is much less than what is necessary when using a 2nd order approach. However, the computational cost is significantly reduced when considering 1st order derivatives or 2nd order derivatives without cross-terms. In the following, we will only consider the empirical, first order and 2nd order without cross-terms approaches.

In Figs 4 and A4, we compare the diagonal of the $C_p$ matrix for HRGPS and strong motion stations. The 1st and 2nd order matrices seem to capture the main features of the empirical $C_p$ matrix. Overall, the diagonal elements of the 2nd order $C_p$ are more similar to the empirical $C_p$ matrix. Even if the 2nd order $C_p$ is computed after neglecting 2nd order cross-terms in eq. (9), Fig. A5 shows that the difference with respect to the empirical $C_p$ matrix is 10–20 per cent smaller than the 1st order $C_p$ matrix. Such differences could impact the inversion results. For this reason, in the next section, we explore the impact of the type of $C_p$ matrix estimate on the coseismic models of the 2015 $M_w = 8.3$ Illapel earthquake.

To model the 2015 $M_w = 8.3$ Illapel earthquake, we design a curved fault geometry using the GOCAD® commercial software package matching local seismicity and aftershock focal mechanisms (Fig. 1). The focal mechanisms are from Global CMT (Dziewonski *et al.* 1981) over a period of one month after the main shock. The fault surface is divided into 10 patches along-dip and 17 patches along-strike (170 in total) with 18 km side-length, which in a sense,
<figure><img src="image_6.png" /><figcaption>**Figure 3.** Model variability of the *P* wave, *S* wave and density as a function of depth in the Illapel region. The black line represents the velocity layered model used for Green's function (GF) calculation. Grey histograms are the probability density function for each parameter as a function of depth.</figcaption></figure>

<table><caption>**Table 1.** Approaches to calculate $C_p$ (for 36 parameters).</caption><thead><tr><th>Approach</th><th>Number of forward model evaluations</th></tr></thead><tbody><tr><td>Without C<sub>p</sub></td><td>0</td></tr><tr><td>Empirical</td><td>195 (in this study)</td></tr><tr><td>1st order forward derivatives</td><td>37</td></tr><tr><td>1st order centred derivatives</td><td>72</td></tr><tr><td>2nd order without cross-terms</td><td>73</td></tr><tr><td>2nd order</td><td>1333</td></tr></tbody></table>

We follow three different schemes to map earth model uncertainty into prediction uncertainty. The first straightforward approach is to empirically calculate the prediction uncertainty covariance matrix $C_p$ using predictions computed for a large number of random earth models $\Psi^i$, ($i = 1, \dots, n$) drawn from $p(\log \Psi)$:

$$C_p = \frac{1}{n-1} \sum_{i=1}^{n} (g(\Psi^i, \mathbf{m}) - g(\tilde{\Psi}, \mathbf{m})) (g(\Psi^i, \mathbf{m}) - g(\tilde{\Psi}, \mathbf{m}))^T, (5)$$

where $g(\Psi^i, \mathbf{m})$ is the prediction for the earth model $\Psi^i$ and the source model $\mathbf{m}$. In our case, we use a preliminary source model $\mathbf{m}$ derived from a first preliminary slip inversion. $g(\tilde{\Psi}, \mathbf{m})$ is the prediction response for the average earth model $\tilde{\Psi}$. This empirical approach is computationally expensive because it needs the calculation of predictions for each randomly generated earth model. To evaluate the number of models $n$ necessary to calculate an accurate empirical $C_p$ matrix, we compare empirical $C_p$ matrices calculated for an increasing number of random earth models. We observe that the empirical $C_p$ matrix is converging using 195 random Earth samples (Fig. A3), corresponding to relatively smooth histograms in Fig. 3.

To test a computationally less expensive approach, we also follow the first-order approximation approach proposed by Duputel *et al.* (2014). Assuming that we can approximate our forward model $g(\Psi, \mathbf{m})$ by linearized perturbations, for an *a priori* earth model we have

then:

$$g(\Psi, \mathbf{m}) \approx g(\tilde{\Psi}, \mathbf{m}) + K_{\Psi}(\tilde{\Psi}, \mathbf{m}) \cdot (\Psi - \tilde{\Psi}), \quad (6)$$

where $K$ is the sensitivity kernels of the prediction with respect to elastic parameters used to compute forward predictions:

$$(K_{\Psi})_{ij}(\tilde{\Psi}, \mathbf{m}) = \frac{\partial g_i}{\partial \Psi_j}(\tilde{\Psi}, \mathbf{m}), \quad (7)$$

where $\Psi_j$ corresponds to the jth elastic parameter in the earth model $\Psi$. We use then $K$ to estimate $C_p$ as:

$$C_p = K_{\Psi} \cdot C_{\Psi} \cdot K_{\Psi}^{\mathrm{T}}, \quad (8)$$

where $C_{\Psi}$ is the same log-normal covariance that we use for perturbing the random models of the empirical $C_p$ in eq. (4). Although this approach looks appropriate for static data, it could be problematic for kinematic data as the link between earth model perturbations and waveform predictions is probably not linear. Indeed, changes in the velocity model induce both time-shifts and amplitude variations in the predicted waveforms.

Therefore, we also explore the possibility of using a 2nd order perturbation approach of the forward model as:

$$g(\Psi, \mathbf{m}) \approx g(\tilde{\Psi}, \mathbf{m}) + K_{\Psi}(\tilde{\Psi}, \mathbf{m}) \cdot (\Psi - \tilde{\Psi}) + \frac{1}{2} (\Psi - \tilde{\Psi}) \cdot H_{\Psi}(\tilde{\Psi}, \mathbf{m}) \cdot (\Psi - \tilde{\Psi}), \quad (9)$$

where $H_{\Psi}$ includes the second order derivative with respect to the elastic parameters:

$$(H_{\Psi})_{ijk}(\tilde{\Psi}, \mathbf{m}) = \frac{\partial^2 g_i}{\partial \Psi_k \partial \Psi_j}(\tilde{\Psi}, \mathbf{m}). \quad (10)$$

From eq. (9), we can then calculate the $C_p$ matrix using eq. (5) by rapidly generating a large number of forward model predictions.

The derivatives in eq. (9) are computed numerically using finite differences. We summarize the difference in computational cost between approaches in table 1. The computational cost of each approach in terms of forward model evaluation is summarized in Table 1. In this study, the empirical approach necessitated about 200 forward model evaluations, which is much less than what is necessary when using a 2nd order approach. However, the computational cost is significantly reduced when considering 1st order derivatives or 2nd order derivatives without cross-terms. In the following, we will only consider the empirical, first order and 2nd order without cross-terms approaches.

In Figs 4 and A4, we compare the diagonal of the $C_p$ matrix for HRGPS and strong motion stations. The 1st and 2nd order matrices seem to capture the main features of the empirical $C_p$ matrix. Overall, the diagonal elements of the 2nd order $C_p$ are more similar to the empirical $C_p$ matrix. Even if the 2nd order $C_p$ is computed after neglecting 2nd order cross-terms in eq. (9), Fig. A5 shows that the difference with respect to the empirical $C_p$ matrix is 10–20 per cent smaller than the 1st order $C_p$ matrix. Such differences could impact the inversion results. For this reason, in the next section, we explore the impact of the type of $C_p$ matrix estimate on the coseismic models of the 2015 $M_w = 8.3$ Illapel earthquake.

To model the 2015 $M_w = 8.3$ Illapel earthquake, we design a curved fault geometry using the GOCAD® commercial software package matching local seismicity and aftershock focal mechanisms (Fig. 1). The focal mechanisms are from Global CMT (Dziewonski *et al.* 1981) over a period of one month after the main shock. The fault surface is divided into 10 patches along-dip and 17 patches along-strike (170 in total) with 18 km side-length, which in a sense,
<figure><img src="image_7.png" /><figcaption>(a) HRGPS</figcaption></figure>

<figure><img src="image_8.png" /><figcaption>(b) Strong motion</figcaption></figure>

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

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

**Figure 4.** Covariance matrix comparison for HRGPS records (a) and strong motion stations (b). The green line represents the diagonal of the empirical covariance matrix (i.e. the matrix created from an ensemble of models). The red and blue line represents the diagonal of the matrix calculated using the 1st and 2nd order approximation approach.

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

**Figure 5.** Posterior mean coseismic slip model for the static data set. Arrows represent the slip directions and the ellipses their associated uncertainties assuming a 95 per cent confidence interval.

is a spatial regularization. However, we do not impose any smoothing or empirical regularizations in the inverse problem, which could potentially smooth out rupture complexities. For the static inversion, we invert for along-strike and along-dip slip components in each subfault. In the full joint inversion, we invert for both slip components along with rise time, rupture velocity, and the hypocentre location on the fault (along-strike and along-dip distance). We model the rupture front by solving the eikonal equation for a candidate rupture velocity in each subfault. Each subfault is discretized into 10 × 10 point sources that rupture sequentially as the rupture front passes. During the earthquake, each point on the fault is allowed to rupture only once [contrary to a multiwindow approach; Hartzell & Heaton (1983), Li *et al.* (2016)], adopting a prescribed triangular slip rate function. Even though multi-window approach is able to recover great complexity in the slip rate functions, the single window approach works better for recovering rupture velocity and seismic moment and at the same time, it significantly decreases the number of inverted parameters (Cohee & Beroza 1994).

In the Bayesian inversion approach, we describe a *priori* PDFs to represent our prior knowledge for each of the parameters to invert. The corresponding *a priori* distributions of our joint model are shown in Fig. A6. We use the hypocentre of the CSN as *a priori* since it was obtained using regional data. For InSAR images, we include a nuisance parameter to correct each image from a constant offset (i.e. two nuisance parameters in total), and for the GPS data sets we add translation parameters (i.e. three parameters for each
<figure><img src="image_7.png" /><figcaption>(a) HRGPS</figcaption></figure>

<figure><img src="image_8.png" /><figcaption>(b) Strong motion</figcaption></figure>

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

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

**Figure 4.** Covariance matrix comparison for HRGPS records (a) and strong motion stations (b). The green line represents the diagonal of the empirical covariance matrix (i.e. the matrix created from an ensemble of models). The red and blue line represents the diagonal of the matrix calculated using the 1st and 2nd order approximation approach.

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

**Figure 5.** Posterior mean coseismic slip model for the static data set. Arrows represent the slip directions and the ellipses their associated uncertainties assuming a 95 per cent confidence interval.

is a spatial regularization. However, we do not impose any smoothing or empirical regularizations in the inverse problem, which could potentially smooth out rupture complexities. For the static inversion, we invert for along-strike and along-dip slip components in each subfault. In the full joint inversion, we invert for both slip components along with rise time, rupture velocity, and the hypocentre location on the fault (along-strike and along-dip distance). We model the rupture front by solving the eikonal equation for a candidate rupture velocity in each subfault. Each subfault is discretized into 10 × 10 point sources that rupture sequentially as the rupture front passes. During the earthquake, each point on the fault is allowed to rupture only once [contrary to a multiwindow approach; Hartzell & Heaton (1983), Li *et al.* (2016)], adopting a prescribed triangular slip rate function. Even though multi-window approach is able to recover great complexity in the slip rate functions, the single window approach works better for recovering rupture velocity and seismic moment and at the same time, it significantly decreases the number of inverted parameters (Cohee & Beroza 1994).

In the Bayesian inversion approach, we describe *a priori* PDFs to represent our prior knowledge for each of the parameters to invert. The corresponding *a priori* distributions of our joint model are shown in Fig. A6. We use the hypocentre of the CSN as *a priori* since it was obtained using regional data. For InSAR images, we include a nuisance parameter to correct each image from a constant offset (i.e. two nuisance parameters in total), and for the GPS data sets we add translation parameters (i.e. three parameters for each
<figure><img src="image_12.png" /><figcaption>**Figure 6.** Comparison of coseismic slip distributions obtained using different prediction error covariances $C_p$. Red colours indicate the corresponding posterior mean coseismic slip model. Arrows represent the slip directions with their corresponding uncertainty shown by ellipses. The red star is the inverted hypocentre location (empirical, 1st and 2nd order approximation, respectively). The blue star is the CSN hypocentre, and the green star is the USGS hypocentre.</figcaption></figure>

<figure><img src="image_13.png" /><figcaption>**Figure 7.** Rupture times comparison between different $C_p$ inversion solutions. Comparison between the empirical covariance matrix and the first order (a) and 2nd order (b) approaches.</figcaption></figure>

set). These parameters are used to redefine the reference frame of each geodetic data set during the inversion process, since both InSAR and GNSS are relative measurements, and have their own reference frame.

Since we are working with different data sets, we want to know how sensitive they are to slip on the fault. Thus, we carry out a sensitivity analysis for each data set. We follow an approach similar to Duputel *et al.* (2015). The sensitivity of each data set is calculated as:

$$S(D) = \operatorname{diag}(G'(D) \cdot C_{\chi}^{-1}(D) \cdot G(D)), \qquad (11)$$

where $G$ is the corresponding Green functions (in the along-dip direction), and $C_\chi$ is the covariance matrix described above for a given data set $D$. For a given subfault, this measure is equivalent to computing the $L_2$ norm of the predictions due to unit dip-slip in the considered patch. The corresponding sensitivities are shown in Fig. A7. GPS and InSAR data sets are sensitive to slip in most fault areas, except for the shallowest region. On the other hand, tsunami data is not sensitive to slip in the inshore fault region but to the offshore zone. The kinematic data is globally sensitive to slip over the entire fault. Finally, if we use the whole data set, although we still observe a decrease in sensitivity at the
<figure><img src="image_12.png" /><figcaption>**Figure 6.** Comparison of coseismic slip distributions obtained using different prediction error covariances $C_p$. Red colours indicate the corresponding posterior mean coseismic slip model. Arrows represent the slip directions with their corresponding uncertainty shown by ellipses. The red star is the inverted hypocentre location (empirical, 1st and 2nd order approximation, respectively). The blue star is the CSN hypocentre, and the green star is the USGS hypocentre.</figcaption></figure>

<figure><img src="image_13.png" /><figcaption>**Figure 7.** Rupture times comparison between different $C_p$ inversion solutions. Comparison between the empirical covariance matrix and the first order (a) and 2nd order (b) approaches.</figcaption></figure>

set). These parameters are used to redefine the reference frame of each geodetic data set during the inversion process, since both InSAR and GNSS are relative measurements, and have their own reference frame.

Since we are working with different data sets, we want to know how sensitive they are to slip on the fault. Thus, we carry out a sensitivity analysis for each data set. We follow an approach similar to Duputel et al. (2015). The sensitivity of each data set is calculated as:

$$S(D) = \operatorname{diag}(G'(D) \cdot C_{\chi}^{-1}(D) \cdot G(D)), \qquad (11)$$

where G is the corresponding Green functions (in the along-dip direction), and $C_\chi$ is the covariance matrix described above for a given data set D. For a given subfault, this measure is equivalent to computing the $L_2$ norm of the predictions due to unit dip-slip in the considered patch. The corresponding sensitivities are shown in Fig. A7. GPS and InSAR data sets are sensitive to slip in most fault areas, except for the shallowest region. On the other hand, tsunami data is not sensitive to slip in the inshore fault region but to the offshore zone. The kinematic data is globally sensitive to slip over the entire fault. Finally, if we use the whole data set, although we still observe a decrease in sensitivity at the
<figure><img src="image_14.png" /><figcaption>(b) Moment Rate functions</figcaption></figure>

<figure><img src="image_15.png" /><figcaption>(c) Centroid times</figcaption></figure>

<figure><img src="image_16.png" /><figcaption>(d) Slip Uncertainty</figcaption></figure>

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

**Figure 8.** Impact of using a 2nd order approximation $C_p$ in slip inversion. (a) Posterior mean coseismic slip model, arrows represent the slip directions and the ellipses its corresponding uncertainty. Contours show stochastic rupture fronts samples from the *a posteriori* distribution every 10 s. (b) Stochastic moment rate functions. (c) Posterior mean coseismic slip model with contours that represent stochastic centroid time fronts samples from the *a posteriori* distribution. (d) Uncertainty of the ensemble of coseismic slip models. The red star in the figures represents the inverted hypocentre location.

trench, we have an overall good sensitivity to slip over the entire fault.

## 4 RESULTS

According to our cascading approach, we first perform an inversion of the final slip using static data (that is, InSAR, GPS and tsunami data). We thus generate a posterior ensemble of slip models whose posterior mean and uncertainty is shown in Fig. 5. This model presents two main slip patches that extend updip to the trench. The solution obtained using static data only has a peak slip of about 10.9 ± 16.0 m, while the mean fault slip is about 2.5 ± 1.8 m (assuming a 95 per cent confidence interval). We observe that uncertainties are as large as the posterior mean slip amplitude. In addition, we see that even if tsunami data is employed, slip uncertainty is larger in the shallow part of the fault, due to the lack of data coverage in that area.

We then use the *a posteriori* PDF of the static slip model as a starting point to make three different joint inversions: (i) a joint

inversion using an empirical $C_p$ matrix, (ii) a joint inversion using a $C_p$ matrix calculated using the first-order perturbation approach and (iii) a joint inversion using a $C_p$ matrix calculated using the second-order perturbation approach. The posterior mean coseismic slip models obtained using these different approaches are shown in Fig. 6. We also compare the posterior distributions of dip-slip in the online supplement (Fig. A8). The three solutions exhibit two principal slip regions, one northwestward of the hypocentre and another at shallow depth reaching the trench. The deeper slip patch is well constrained for the three solutions, with a mean slip of 6 to meters for this region. The solution based on 1st order $C_p$ shows a compact slip patch at shallow depth, while shallow slip is more broadly distributed when considering 2nd order or empirical $C_p$ matrices. This results into a larger peak slip value for the 1st order $C_p$ solution (21.0 ± 4.1 m), while solutions obtained with an empirical $C_p$ (15.88 ± 5.0 m) and with a 2nd order $C_p$ (17.63 ± 6.8 m) display smaller peak slip values. Uncertainties significantly decrease when incorporating the kinematic data set.
<figure><img src="image_14.png" /><figcaption>(b) Moment Rate functions</figcaption></figure>

<figure><img src="image_15.png" /><figcaption>(c) Centroid times</figcaption></figure>

<figure><img src="image_16.png" /><figcaption>(d) Slip Uncertainty</figcaption></figure>

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

**Figure 8.** Impact of using a 2nd order approximation $C_p$ in slip inversion. (a) Posterior mean coseismic slip model, arrows represent the slip directions and the ellipses its corresponding uncertainty. Contours show stochastic rupture fronts samples from the *a posteriori* distribution every 10 s. (b) Stochastic moment rate functions. (c) Posterior mean coseismic slip model with contours that represent stochastic centroid time fronts samples from the *a posteriori* distribution. (d) Uncertainty of the ensemble of coseismic slip models. The red star in the figures represents the inverted hypocentre location.

trench, we have an overall good sensitivity to slip over the entire fault.

## 4 RESULTS

According to our cascading approach, we first perform an inversion of the final slip using static data (that is, InSAR, GPS and tsunami data). We thus generate a posterior ensemble of slip models whose posterior mean and uncertainty is shown in Fig. 5. This model presents two main slip patches that extend updip to the trench. The solution obtained using static data only has a peak slip of about 10.9 ± 16.0 m, while the mean fault slip is about 2.5 ± 1.8 m (assuming a 95 per cent confidence interval). We observe that uncertainties are as large as the posterior mean slip amplitude. In addition, we see that even if tsunami data is employed, slip uncertainty is larger in the shallow part of the fault, due to the lack of data coverage in that area.

We then use the *a posteriori* PDF of the static slip model as a starting point to make three different joint inversions: (i) a joint

inversion using an empirical $C_p$ matrix, (ii) a joint inversion using a $C_p$ matrix calculated using the first-order perturbation approach and (iii) a joint inversion using a $C_p$ matrix calculated using the second-order perturbation approach. The posterior mean coseismic slip models obtained using these different approaches are shown in Fig. 6. We also compare the posterior distributions of dip-slip in the online supplement (Fig. A8). The three solutions exhibit two principal slip regions, one northwestward of the hypocentre and another at shallow depth reaching the trench. The deeper slip patch is well constrained for the three solutions, with a mean slip of 6 to meters for this region. The solution based on 1st order $C_p$ shows a compact slip patch at shallow depth, while shallow slip is more broadly distributed when considering 2nd order or empirical $C_p$ matrices. This results into a larger peak slip value for the 1st order $C_p$ solution (21.0 ± 4.1 m), while solutions obtained with an empirical $C_p$ (15.88 ± 5.0 m) and with a 2nd order $C_p$ (17.63 ± 6.8 m) display smaller peak slip values. Uncertainties significantly decrease when incorporating the kinematic data set.
<figure><img src="image_18.png" /><figcaption>**Figure 9.** (a) Observed horizontal GPS (black arrows) and predictions for the posterior mean model (red arrows) using a 2nd order approximation $C_p$. (b) The colourmap indicates vertical component displacements for observed GPS (outer circle) and vertical predictions for the posterior mean model (inner circle).</figcaption></figure>

Fig. 7 compares rupture times between solutions (taking the solution based on empirical $C_p$ as reference). Both models obtained using a first and second order $C_p$ result in rupture times similar to those obtained with an empirical covariance matrix. However, the second order approach presents an overall smaller dispersion ($\sigma = 4.75$ s) compared to the first order approach ($\sigma = 5.97$ s). Regardless of the prediction error covariance matrix, we note that the nuisance parameters associated with GPS data sets converge to zero, which means they don't need further corrections (Fig. A9). There is no significant variation in the constant offset associated with the descending InSAR image, with a posterior mean value of 3.7 cm. On the other hand, there are some differences in the nuisance parameter of the ascending interferogram, which vary between -2.5 and -1.5 cm between the different solutions.

Details of the solution obtained using a 2nd order $C_p$ are shown in Fig. 8. Similar figures are presented for the 1st order and empirical $C_p$ in supplementary Figs A10 and A11, respectively. Stochastic rupture propagation fronts in Fig. 8(a) suggest a complex rupture pattern. It slowly grows close to the hypocentre, and then propagates updip, with a rupture speed from 2 to 4 km s⁻¹. Stochastic moment rate functions in Fig. 8(b) indicate an overall rupture duration of 120 s approximately. The average scalar seismic moment is $M_0 = 3.20 \pm 0.045 \times 10^{21}$ Nm, that is a moment magnitude of $M_w = 8.27 \pm 0.005$. We can notice two energy peaks, a small one at 25 s, and another one at 50 s. As it has been reported before (Gombert et al. 2018b), we observe a negative correlation of rise time and initial rupture times (Fig. A12a). However, this correlation disappears when comparing rise time and slip pulse centroid times (Fig. A12b). This arises from the fact that observations are more sensitive to the slip pulse centroid time at each subfault, rather than the initial rupture time and rise time (see Fig. 7; Gombert et al. 2018b). The distribution of centroid times in Fig. 8(c) shows a heterogeneous rupture propagation. In particular, there are regions at

the northwest of the hypocentre that break faster than their corresponding adjacent areas. These complexities are discussed further in Section 5.2.

We use the posterior coseismic model to calculate synthetic displacements and compare them to GPS observations (Fig. 9). Both permanent stations and campaign survey stations show an acceptable fit, including the vertical components. The corresponding residuals are shown in Fig. A13. The residuals are globally small compared with uncertainties. For the horizontal components, the average residual is approximately 10 cm, while for the vertical component is 5 cm, which is acceptable given the magnitude of the displacements (up to 2 m). Stochastic predictions of tsunami waveforms display a good agreement with tsunami observations (Fig. 10). In particular, we see that later arrivals are often well fitted even if they are not included in the data set used for the slip inversion. The tide gauges buca1, papol and talt1, and the DART stations D32411, D43412 and D51407 present a slight time-shift between observed and predicted waveforms. This shift could be explained by local site effects, local bathymetry for the case of tide gauges, and in the case of DART stations, by path trajectory not accurately modelled by the forward model. Fig. 11 shows that InSAR data is also well predicted by our posterior coseismic model, with residuals smaller than 10 per cent of maximum LOS displacements. The spatial distribution of the residuals does not correlate with the coseismic displacement pattern. Nevertheless, we notice a spatial pattern in the InSAR residuals from Fig. 11(c), with more positive values in the northern region of the image. This residual could be linked to discrepancies between different types of geodetic observations in the region. Although the ascending image shows vertical displacements up to 40 cm, GPS vertical displacements in the same area are close to zero or even display negative values. This InSAR residual pattern has also been observed by Klein et al. (2017). We also use the posterior coseismic model to calculate kinematic stochastic waveforms. Kinematic data show a directivity effect with larger amplitudes towards the north
<figure><img src="image_18.png" /><figcaption>**Figure 9.** (a) Observed horizontal GPS (black arrows) and predictions for the posterior mean model (red arrows) using a 2nd order approximation $C_p$. (b) The colourmap indicates vertical component displacements for observed GPS (outer circle) and vertical predictions for the posterior mean model (inner circle).</figcaption></figure>

Fig. 7 compares rupture times between solutions (taking the solution based on empirical $C_p$ as reference). Both models obtained using a first and second order $C_p$ result in rupture times similar to those obtained with an empirical covariance matrix. However, the second order approach presents an overall smaller dispersion ($\sigma = 4.75$ s) compared to the first order approach ($\sigma = 5.97$ s). Regardless of the prediction error covariance matrix, we note that the nuisance parameters associated with GPS data sets converge to zero, which means they don't need further corrections (Fig. A9). There is no significant variation in the constant offset associated with the descending InSAR image, with a posterior mean value of 3.7 cm. On the other hand, there are some differences in the nuisance parameter of the ascending interferogram, which vary between -2.5 and -1.5 cm between the different solutions.

Details of the solution obtained using a 2nd order $C_p$ are shown in Fig. 8. Similar figures are presented for the 1st order and empirical $C_p$ in supplementary Figs A10 and A11, respectively. Stochastic rupture propagation fronts in Fig. 8(a) suggest a complex rupture pattern. It slowly grows close to the hypocentre, and then propagates updip, with a rupture speed from 2 to 4 km s⁻¹. Stochastic moment rate functions in Fig. 8(b) indicate an overall rupture duration of 120 s approximately. The average scalar seismic moment is $M_0 = 3.20 \pm 0.045 \times 10^{21}$ Nm, that is a moment magnitude of $M_w = 8.27 \pm 0.005$. We can notice two energy peaks, a small one at 25 s, and another one at 50 s. As it has been reported before (Gombert et al. 2018b), we observe a negative correlation of rise time and initial rupture times (Fig. A12a). However, this correlation disappears when comparing rise time and slip pulse centroid times (Fig. A12b). This arises from the fact that observations are more sensitive to the slip pulse centroid time at each subfault, rather than the initial rupture time and rise time (see Fig. 7; Gombert et al. 2018b). The distribution of centroid times in Fig. 8(c) shows a heterogeneous rupture propagation. In particular, there are regions at

the northwest of the hypocentre that break faster than their corresponding adjacent areas. These complexities are discussed further in Section 5.2.

We use the posterior coseismic model to calculate synthetic displacements and compare them to GPS observations (Fig. 9). Both permanent stations and campaign survey stations show an acceptable fit, including the vertical components. The corresponding residuals are shown in Fig. A13. The residuals are globally small compared with uncertainties. For the horizontal components, the average residual is approximately 10 cm, while for the vertical component is 5 cm, which is acceptable given the magnitude of the displacements (up to 2 m). Stochastic predictions of tsunami waveforms display a good agreement with tsunami observations (Fig. 10). In particular, we see that later arrivals are often well fitted even if they are not included in the data set used for the slip inversion. The tide gauges buca1, papol and talt1, and the DART stations D32411, D43412 and D51407 present a slight time-shift between observed and predicted waveforms. This shift could be explained by local site effects, local bathymetry for the case of tide gauges, and in the case of DART stations, by path trajectory not accurately modelled by the forward model. Fig. 11 shows that InSAR data is also well predicted by our posterior coseismic model, with residuals smaller than 10 per cent of maximum LOS displacements. The spatial distribution of the residuals does not correlate with the coseismic displacement pattern. Nevertheless, we notice a spatial pattern in the InSAR residuals from Fig. 11(c), with more positive values in the northern region of the image. This residual could be linked to discrepancies between different types of geodetic observations in the region. Although the ascending image shows vertical displacements up to 40 cm, GPS vertical displacements in the same area are close to zero or even display negative values. This InSAR residual pattern has also been observed by Klein et al. (2017). We also use the posterior coseismic model to calculate kinematic stochastic waveforms. Kinematic data show a directivity effect with larger amplitudes towards the north
<figure><img src="image_19.png" /><figcaption>**Figure 10.** Comparisons between tsunami observations (black) and stochastic predictions (red) using a 2nd order approximation $C_p$. The tsunami waveform signal used in the inversion is shown between blue dots. The map depicts each tsunami station locations.</figcaption></figure>

that is well reproduced by the model (Figs 12 and A14). We can see that stochastic waveforms reproduce most of the features visible in the HRGPS and strong motion records, even at large distances (i.e. distances >2°).

## 5 DISCUSSION

We compare our slip models with previous models published in the literature. Our posterior coseismic model presents a maximum slip of 17.63 ± 6.8 m at shallow depth. This slip magnitude is larger than the one observed by Klein et al. (2017) (10 m); Ruiz et al. (2016), Shrivastava et al. (2016), An & Meng (2017) (8 ms) and previous kinematic models such as the one of Tilmann et al. (2016), Heidarzadeh et al. (2016), Li et al. (2016), Melgar et al. (2016) (6–12 m). Overall, our joint model is more similar to the slip distribution of Melgar et al. (2016), which exhibits two slip regions, with a maximum slip of 12 m, which is smaller than our posterior mean estimate but within uncertainty of our solution. This difference likely results from the fact that our results rely only on spatial discretization in square subfaults while the inversion does not incorporate smoothing constraints, contrary to the aforementioned studies that incorporate smoothing regularizations. By using such constraints, the slip distributions are smoother, which penalizes abrupt changes and locally high slip amplitudes.

Although largest slip amplitudes in our posterior model are located at relatively shallow depth, we note that several previously published models include slip extending to deeper regions of the fault (i.e. below the coast). In this regard, Klein et al. (2017) suggest that slip at larger depth is necessary to fit vertical GPS observations. Although the fault slips mostly offshore according to our solution, we still observe significant slip (2–3 m) at larger depth. In Fig. A15, we investigate the contribution of slip at different depths to fit the vertical GPS observations. In agreement with our sensitivity maps in Fig. A7, we see that shallow slip does not generate much displacement inland. Although we see that a moderate amount of slip close to the coast generates uplift in our model predictions, our model still features some misfits on coastal GPS stations (as shown in Fig. A13), which can explain the difference in the amount of slip at depth compared to previous models (e.g. at station EMAT with an observed uplift of 20 cm, Klein et al. (2017) has a misfit of 5 cm, while our solution corresponds to a misfit of 8 cm).

In the next subsections, we will examine individually different aspects of the Illapel earthquake rupture. We first assess the reliability of our model close to the trench by exploring the importance of shallow slip to fit tsunami records. We then investigate encircling rupture patterns visible in our solutions.
<figure><img src="image_19.png" /><figcaption>**Figure 10.** Comparisons between tsunami observations (black) and stochastic predictions (red) using a 2nd order approximation $C_p$. The tsunami waveform signal used in the inversion is shown between blue dots. The map depicts each tsunami station locations.</figcaption></figure>

that is well reproduced by the model (Figs 12 and A14). We can see that stochastic waveforms reproduce most of the features visible in the HRGPS and strong motion records, even at large distances (i.e. distances >2°).

## 5 DISCUSSION

We compare our slip models with previous models published in the literature. Our posterior coseismic model presents a maximum slip of 17.63 ± 6.8 m at shallow depth. This slip magnitude is larger than the one observed by Klein et al. (2017) (10 m); Ruiz et al. (2016), Shrivastava et al. (2016), An & Meng (2017) (8 ms) and previous kinematic models such as the one of Tilmann et al. (2016), Heidarzadeh et al. (2016), Li et al. (2016), Melgar et al. (2016) (6–12 m). Overall, our joint model is more similar to the slip distribution of Melgar et al. (2016), which exhibits two slip regions, with a maximum slip of 12 m, which is smaller than our posterior mean estimate but within uncertainty of our solution. This difference likely results from the fact that our results rely only on spatial discretization in square subfaults while the inversion does not incorporate smoothing constraints, contrary to the aforementioned studies that incorporate smoothing regularizations. By using such constraints, the slip distributions are smoother, which penalizes abrupt changes and locally high slip amplitudes.

Although largest slip amplitudes in our posterior model are located at relatively shallow depth, we note that several previously published models include slip extending to deeper regions of the fault (i.e. below the coast). In this regard, Klein et al. (2017) suggest that slip at larger depth is necessary to fit vertical GPS observations. Although the fault slips mostly offshore according to our solution, we still observe significant slip (2–3 m) at larger depth. In Fig. A15, we investigate the contribution of slip at different depths to fit the vertical GPS observations. In agreement with our sensitivity maps in Fig. A7, we see that shallow slip does not generate much displacement inland. Although we see that a moderate amount of slip close to the coast generates uplift in our model predictions, our model still features some misfits on coastal GPS stations (as shown in Fig. A13), which can explain the difference in the amount of slip at depth compared to previous models (e.g. at station EMAT with an observed uplift of 20 cm, Klein et al. (2017) has a misfit of 5 cm, while our solution corresponds to a misfit of 8 cm).

In the next subsections, we will examine individually different aspects of the Illapel earthquake rupture. We first assess the reliability of our model close to the trench by exploring the importance of shallow slip to fit tsunami records. We then investigate encircling rupture patterns visible in our solutions.
<figure><img src="image_20.png" /><figcaption>**Figure 11.** InSAR misfit using the posterior coseismic model using the 2nd order C<sub>p</sub> matrix solution. Observed ascending (a) and descending (d) Sentinel-1A images. We show the corresponding synthetic displacement for ascending (b) and descending (c) images and the respective residual, (c) for ascending and (f) for descending images. Small black arrows represent the LOS and orbit direction, respectively.</figcaption></figure>

## 5.1 Impact of shallow slip

At present, there is no general agreement regarding the amount of shallow slip during the Illapel earthquake since some studies indicate the absence of shallow slip (An & Meng 2017), while others demonstrate that shallow slip is necessary to explain tsunami observations (Lay et al. 2016). To analyse the amount of shallow slip, we evaluate the cumulative posterior PDF of slip in the shallow region (Fig. 13). We observe that the probability of slip to be greater than 13 m at shallow depth is about 83.8 per cent.

To further explore the contribution of shallow slip, we perform a static slip inversion imposing shallow slip to be very small (i.e. in the two shallowest subfault rows). The aforementioned was performed by fixing a prior PDF with a narrow Gaussian centred on zero for the along-dip component of slip (considering a standard deviation of 0.5 m). The corresponding posterior mean model is shown in Fig. 14. If we compare the resulting solution in Fig. 14 with the previous posterior coseismic models in Figs 5 and 8, we can still

find the slip patch close to the hypocentre (longitude -72°, latitude -31.25°). However, the shallow part of the model is significantly different due to the new prior. Regarding the data fit, we can notice that GPS fits remain unchanged between static models (Fig. A16; i.e. GPS observations are insensitive to shallow slip). The comparison of model performance for tsunami observations for both solutions is shown in Fig. 14(b). We notice that the RMS misfit for tsunami data are smaller when including shallow slip (Fig. A17). However, such comparison can be misleading: the model with shallow slip will naturally better fit the observations as it includes more free parameters than the one for which shallow slip is proscribed. To evaluate if the decrease in tsunami misfit is significant, we evaluate two different information criteria: The Bayesian Information Criterion (BIC) and the Akaike Information Criterion (AIC; Bishop 2006) (Supporting Information text S5). In Table 2 we show the differences ΔBIC and ΔAIC, with respect to our solution including shallow slip. Both criteria tend to favor occurrence of shallow slip
<figure><img src="image_20.png" /><figcaption>**Figure 11.** InSAR misfit using the posterior coseismic model using the 2nd order C<sub>p</sub> matrix solution. Observed ascending (a) and descending (d) Sentinel-1A images. We show the corresponding synthetic displacement for ascending (b) and descending (c) images and the respective residual, (c) for ascending and (f) for descending images. Small black arrows represent the LOS and orbit direction, respectively.</figcaption></figure>

## 5.1 Impact of shallow slip

At present, there is no general agreement regarding the amount of shallow slip during the Illapel earthquake since some studies indicate the absence of shallow slip (An & Meng 2017), while others demonstrate that shallow slip is necessary to explain tsunami observations (Lay et al. 2016). To analyse the amount of shallow slip, we evaluate the cumulative posterior PDF of slip in the shallow region (Fig. 13). We observe that the probability of slip to be greater than 13 m at shallow depth is about 83.8 per cent.

To further explore the contribution of shallow slip, we perform a static slip inversion imposing shallow slip to be very small (i.e. in the two shallowest subfault rows). The aforementioned was performed by fixing a prior PDF with a narrow Gaussian centred on zero for the along-dip component of slip (considering a standard deviation of 0.5 m). The corresponding posterior mean model is shown in Fig. 14. If we compare the resulting solution in Fig. 14 with the previous posterior coseismic models in Figs 5 and 8, we can still

find the slip patch close to the hypocentre (longitude –72°, latitude –31.25°). However, the shallow part of the model is significantly different due to the new prior. Regarding the data fit, we can notice that GPS fits remain unchanged between static models (Fig. A16; i.e. GPS observations are insensitive to shallow slip). The comparison of model performance for tsunami observations for both solutions is shown in Fig. 14(b). We notice that the RMS misfit for tsunami data are smaller when including shallow slip (Fig. A17). However, such comparison can be misleading: the model with shallow slip will naturally better fit the observations as it includes more free parameters than the one for which shallow slip is proscribed. To evaluate if the decrease in tsunami misfit is significant, we evaluate two different information criteria: The Bayesian Information Criterion (BIC) and the Akaike Information Criterion (AIC; Bishop 2006) (Supporting Information text S5). In Table 2 we show the differences ΔBIC and ΔAIC, with respect to our solution including shallow slip. Both criteria tend to favor occurrence of shallow slip
<figure><img src="image_21.png" /><figcaption>**Figure 12.** Examples of comparisons between data (black) and stochastic predictions (red) for HRGPS and Strong Motion stations using a 2nd order approximation $C_p$. On the maps, the blue star represents the hypocentre while circles indicate station location (orange for the station depicted and yellow for the other stations). $\phi$ and $\Delta$ represent the azimuth and distance from the epicentre. The angle $\alpha$ is the horizontal component azimuth (0° north, 90° east).</figcaption></figure>

rather than the solution without slip at shallow depth (i.e. the model with shallow slip is associated with smaller BIC and AIC values). In other words, the difference in RMS misfit is sufficient to justify the existence of slip at shallow depth. It is worth mentioning that tsunami data is the only data set controlling the slip at shallow depth since is the most sensitivity data to this feature (as shown in the sensitivities in Fig. A7. The differences with previous back-projection studies come from the fact that such shallow features are difficult to resolve only using seismic information [as pointed out by Lay *et al.* (2016)].

Finally, we compare the posterior mean joint coseismic slip distribution with aftershocks locations (Fig. 15). We observe aftershocks in the outer-trench zone, distributed along the shallow slip region revealed by our solution. As suggested by Sladen & Trevisan (2018), the occurrence of outer-rise aftershocks can be used as a proxy to estimate the occurrence of slip at shallow depth along the subduction interface. The distribution of aftershocks is therefore consistent with the occurrence of shallow slip during the Illapel earthquake.

The existence of large slip at shallow depth supports the fact that the 2015 event is not a simple repeat of the earthquake that affected the region in 1943 (Tilmann et al. 2016). This is consistent with historical reports indicating that the tsunami generated in 1943 was much smaller than what was observed in 2015. In addition, the differences in the duration of teleseismic body-wave arrivals for both events suggest that the 1943 rupture did not involve shallow slip (Tilmann et al. 2016). The reason why the 2015 event involves

shallow slip contrarily to the 1943 event is unclear. One possibility is that shallow slip deficit was larger in 2015 than in 1943. This is consistent with coupling models from Métois et al. (2016) showing that the fault is not creeping at plate rate at shallow depth. However, this remains speculative as fault coupling close to the trench is poorly resolved by land-based geodetic data and could potentially be biased when ignoring stress shadowing effects (Lindsey et al. 2021).

## 5.2 Encircling rupture pattern during the 2015 Illapel earthquake

Back-projection results from Meng et al. (2018) show an encircling rupture during the 2015 Illapel earthquake. However, this encircling effect has not been reported by any previous kinematic slip inversion model. Results in Figs 8(a) and (c) show a possible encircling behaviour northwestward from the hypocentre location. We use the posterior coseismic mean model to investigate the slip and slip rate evolution. Snapshots from the slip rate history (Supplementary Movie 1) and slip history (Supplementary Movie 2) are shown in Figs 16 and A18, respectively. The rupture slowly grows propagating up-dip for 38 s. During this first stage of the rupture, we observe two different slip rate patches in Supplementary Movie 2, a main region in the updip fault area, and at 32 s, a secondary slip rate patch in the downdip region. This secondary patch rapidly vanishes after a few seconds, without producing significant slip.
<figure><img src="image_21.png" /><figcaption>**Figure 12.** Examples of comparisons between data (black) and stochastic predictions (red) for HRGPS and Strong Motion stations using a 2nd order approximation $C_p$. On the maps, the blue star represents the hypocentre while circles indicate station location (orange for the station depicted and yellow for the other stations). $\phi$ and $\Delta$ represent the azimuth and distance from the epicentre. The angle $\alpha$ is the horizontal component azimuth (0° north, 90° east).</figcaption></figure>

rather than the solution without slip at shallow depth (i.e. the model with shallow slip is associated with smaller BIC and AIC values). In other words, the difference in RMS misfit is sufficient to justify the existence of slip at shallow depth. It is worth mentioning that tsunami data is the only data set controlling the slip at shallow depth since is the most sensitivity data to this feature (as shown in the sensitivities in Fig. A7. The differences with previous back-projection studies come from the fact that such shallow features are difficult to resolve only using seismic information [as pointed out by Lay *et al.* (2016)]).

Finally, we compare the posterior mean joint coseismic slip distribution with aftershocks locations (Fig. 15). We observe aftershocks in the outer-trench zone, distributed along the shallow slip region revealed by our solution. As suggested by Sladen & Trevisan (2018), the occurrence of outer-rise aftershocks can be used as a proxy to estimate the occurrence of slip at shallow depth along the subduction interface. The distribution of aftershocks is therefore consistent with the occurrence of shallow slip during the Illapel earthquake.

The existence of large slip at shallow depth supports the fact that the 2015 event is not a simple repeat of the earthquake that affected the region in 1943 (Tilmann et al. 2016). This is consistent with historical reports indicating that the tsunami generated in 1943 was much smaller than what was observed in 2015. In addition, the differences in the duration of teleseismic body-wave arrivals for both events suggest that the 1943 rupture did not involve shallow slip (Tilmann et al. 2016). The reason why the 2015 event involves

shallow slip contrarily to the 1943 event is unclear. One possibility is that shallow slip deficit was larger in 2015 than in 1943. This is consistent with coupling models from Métois et al. (2016) showing that the fault is not creeping at plate rate at shallow depth. However, this remains speculative as fault coupling close to the trench is poorly resolved by land-based geodetic data and could potentially be biased when ignoring stress shadowing effects (Lindsey et al. 2021).

## 5.2 Encircling rupture pattern during the 2015 Illapel earthquake

Back-projection results from Meng et al. (2018) show an encircling rupture during the 2015 Illapel earthquake. However, this encircling effect has not been reported by any previous kinematic slip inversion model. Results in Figs 8(a) and (c) show a possible encircling behaviour northwestward from the hypocentre location. We use the posterior coseismic mean model to investigate the slip and slip rate evolution. Snapshots from the slip rate history (Supplementary Movie 1) and slip history (Supplementary Movie 2) are shown in Figs 16 and A18, respectively. The rupture slowly grows propagating up-dip for 38 s. During this first stage of the rupture, we observe two different slip rate patches in Supplementary Movie 2, a main region in the updip fault area, and at 32 s, a secondary slip rate patch in the downdip region. This secondary patch rapidly vanishes after a few seconds, without producing significant slip.
<figure><img src="image_22.png" /><figcaption>**Figure 13.** Cumulative probability of having a slip greater or equal to a corresponding amplitude for a subfault experiencing large slip at shallow depth. The corresponding subfault is shown in the inset figure on the left. Colours represent the posterior mean coseismic slip model using the 2nd order approximation approach. Arrows and ellipses represent the slip directions and their corresponding uncertainties, respectively.</figcaption></figure>

Different back-projection studies show a downdip HF source, that radiates energy for at least 60 s (An *et al.* 2017; Melgar *et al.* 2016). Even though the downdip slip rate in our model is only active for 30 s, the location of this patch is similar to the aforementioned back-projection sources. This difference in the duration is probably due to our single window parameterization, since a subfault cannot break several times in our model. However, if we compare the moment rate function of the slip model proposed by An *et al.* (2017) for the updip and downdip regions with the results of Fig. 8(b), we see that we have similar moment rate functions.

Around 40 s after origin time, the rupture separates in three pulses depicting a first encircling pattern updip from the hypocentre and then another encircling pattern above the first one, also updip from the hypocentre (Fig. 16 and Supplementary Movie 2). These encircling slip pulses contour fault areas with smaller slip rates. This is illustrated in Fig. 17 showing the posterior mean peak-slip rates for every point on the fault. All rupture branches finally join together generating a large slip-rate pulse around 60 s, continuing towards the north along the trench until the end of the earthquake. To investigate the reliability of these encircling rupture patterns, we examine the variability of model samples drawn from the posterior PDF. This is shown in the Supplementary Movie 3, which shows the variability of subfault peak slip rates for different samples of our solution. We clearly see that the two encircled regions are consistently surrounded by areas of larger slip rates. This suggests that the two encircling patterns are robust features of our solution.

To identify which part of the waveform is related to the encircled region, we calculate theoretical S wave traveltimes before

and after the first encircled region (Fig. A19). Between these arrival times, we identify a very sharp positive pulse on the east components of stations, in both HRGPS and strong motion, at the north of the hypocentre. This observation is quite consistent with simulations provided by Page *et al.* (2005), which showed that such encircled barriers are associated with sharp secondary pulses in the seismograms. This sharp phase is less visible on southern stations, even if a longer period pulse is visible. This difference probably results from directivity effects, which lead to larger and sharper signals at the northern stations compared to the southern stations.

The rupture complexity of the 2015 Illapel earthquake is highlighted by the heterogeneity in rupture times and local centroid times shown in Fig. 8. To further analyse such complex rupture, we examine two sets of stochastic slip rate functions corresponding to two different regions that present significant slip rate at 45 and 60 s (shown in Fig. A20). Both slip rate functions exhibit maximums that reach more than 1.0 m s⁻¹. The stochastic slip rate functions with a maximum at 45 s correspond to a region in the middle of the fault with a rise-time of about 5 s, while the ones with a maximum at 60 s are for a subfault located at shallow depth with a local rupture duration of 25 s. Some samples of the stochastic slip rate functions at shallow depth begin at the same time and even before the onset of slip in the middle of the fault. Although the peak slip-rates is similar in both fault regions, the longer rise-time at shallow depth results into significant slip close to the trench. Despite such long rupture duration at shallow depth, our model indicate that the rupture follows a pulse-like behaviour
<figure><img src="image_22.png" /><figcaption>**Figure 13.** Cumulative probability of having a slip greater or equal to a corresponding amplitude for a subfault experiencing large slip at shallow depth. The corresponding subfault is shown in the inset figure on the left. Colours represent the posterior mean coseismic slip model using the 2nd order approximation approach. Arrows and ellipses represent the slip directions and their corresponding uncertainties, respectively.</figcaption></figure>

Different back-projection studies show a downdip HF source, that radiates energy for at least 60 s (An *et al.* 2017; Melgar *et al.* 2016). Even though the downdip slip rate in our model is only active for 30 s, the location of this patch is similar to the aforementioned back-projection sources. This difference in the duration is probably due to our single window parameterization, since a subfault cannot break several times in our model. However, if we compare the moment rate function of the slip model proposed by An *et al.* (2017) for the updip and downdip regions with the results of Fig. 8(b), we see that we have similar moment rate functions.

Around 40 s after origin time, the rupture separates in three pulses depicting a first encircling pattern updip from the hypocentre and then another encircling pattern above the first one, also updip from the hypocentre (Fig. 16 and Supplementary Movie 2). These encircling slip pulses contour fault areas with smaller slip rates. This is illustrated in Fig. 17 showing the posterior mean peak-slip rates for every point on the fault. All rupture branches finally join together generating a large slip-rate pulse around 60 s, continuing towards the north along the trench until the end of the earthquake. To investigate the reliability of these encircling rupture patterns, we examine the variability of model samples drawn from the posterior PDF. This is shown in the Supplementary Movie 3, which shows the variability of subfault peak slip rates for different samples of our solution. We clearly see that the two encircled regions are consistently surrounded by areas of larger slip rates. This suggests that the two encircling patterns are robust features of our solution.

To identify which part of the waveform is related to the encircled region, we calculate theoretical S wave traveltimes before

and after the first encircled region (Fig. A19). Between these arrival times, we identify a very sharp positive pulse on the east components of stations, in both HRGPS and strong motion, at the north of the hypocentre. This observation is quite consistent with simulations provided by Page *et al.* (2005), which showed that such encircled barriers are associated with sharp secondary pulses in the seismograms. This sharp phase is less visible on southern stations, even if a longer period pulse is visible. This difference probably results from directivity effects, which lead to larger and sharper signals at the northern stations compared to the southern stations.

The rupture complexity of the 2015 Illapel earthquake is highlighted by the heterogeneity in rupture times and local centroid times shown in Fig. 8. To further analyse such complex rupture, we examine two sets of stochastic slip rate functions corresponding to two different regions that present significant slip rate at 45 and 60 s (shown in Fig. A20). Both slip rate functions exhibit maximums that reach more than 1.0 m s⁻¹. The stochastic slip rate functions with a maximum at 45 s correspond to a region in the middle of the fault with a rise-time of about 5 s, while the ones with a maximum at 60 s are for a subfault located at shallow depth with a local rupture duration of 25 s. Some samples of the stochastic slip rate functions at shallow depth begin at the same time and even before the onset of slip in the middle of the fault. Although the peak slip-rates is similar in both fault regions, the longer rise-time at shallow depth results into significant slip close to the trench. Despite such long rupture duration at shallow depth, our model indicate that the rupture follows a pulse-like behaviour
<figure><img src="image_23.png" /><figcaption>**Figure 14.** (a) Posterior mean coseismic slip model for a static inversion with a non shallow slip *a priori*. Arrows represent the slip directions and the ellipses their associated uncertainties. (b) Comparisons between tsunami observations (black) and stochastic predictions with shallow slip (red) and without shallow slip (blue). The tsunami waveform signal used in the inversion is shown between yellow dots. The map shows the depicted tsunami stations in (a).</figcaption></figure>

<table><caption>**Table 2.** BIC and AIC values with and without shallow slip. Bayesian (BIC) and Akaike (AIC) information criteria are defined in supplementary text S5. ΔBIC and ΔAIC are the difference in BIC and AIC values with respect to the slip model including shallow slip. The values suggest that the shallow slip should be included to properly explain the observations.</caption><thead><tr><td>Model</td><td>ΔBIC</td><td>ΔAIC</td></tr></thead><tbody><tr><td>Shallow slip (348 parameters)</td><td>0</td><td>0</td></tr><tr><td>No shallow slip (314 parameters)</td><td>1001</td><td>1096</td></tr></tbody></table>

since the maximum rise-time is around 30 s, which is considerably shorter than the total rupture time (around 100 s, Heaton 1990).

The encircled slip pulses visible in our solution between 30 and 60 s are consistent with previous back-projection results that suggest such complexities in the rupture (e.g. Ruiz *et al.* 2016; Meng *et al.* 2018). Ruiz *et al.* (2016) show an early stage bilateral rupture that later merged and propagated updip. Meng *et al.* (2018) report two episodes of splitting of rupture fronts, occurring both before reaching 60 s [an effect known as ‘double encircling pincer movement’ (Das & Kostrov 1983)]. The first episode reported by Meng *et al.* (2018) is between 15 and 35 s, and the second, around 45 and 60 s. The first encircling is colocated with the static coseismic model of An & Meng (2017). Consequently, Meng *et al.* (2018) suggest that the encircled region is an asperity. However, this static coseismic model could miss rupture features retrieved by our joint inversion that incorporates additional static and kinematic data. As previously pointed out (Ishii *et al.* 2007; Tilmann *et al.* 2016), back-projection sources trace the progression and changes of the rupture but are not proportional to slip. Our solution is more heterogenous, presenting multiple slip areas with both encircling episodes contouring

regions with small slip rates (and moderate slip), generating particularly high slip rates where the rupture focuses in the final stage of the earthquake (see time=60s, in Figure 16). In this sense, our observations suggest rather the contouring of two regions that do not slip during the rupture. Such strong changes in the rupture propagation associated with high slip rates explain the back-projection results of Meng *et al.* (2018). The small slip amplitude inside the contoured regions can be caused by different factors: (i) these areas could correspond to coupled regions (preventing seismic slip to occur), (ii) complexities at the subduction interface (e.g. due to fracture zones or seamounts) could prevent slip to propagate in these areas or (iii) the contoured regions could be far from the rupture (i.e. initial and dynamic stresses smaller than the fault strength). Regarding the coupling at the subduction interface, the model of the region proposed by Vigny *et al.* (2009) and updated by Métois *et al.* (2012) and Métois *et al.* (2016) shows a relatively high coupling coefficient in the Illapel earthquake area, except in the shallowest region, where the coefficient can be as low as 0.2. However, coupling close to the trench is usually poorly constrained by land-based geodetic data. The Illapel earthquake occurred in the Metropolitan segment defined by Métois *et al.* (2016), and is bounded in the north by the La Serena Low-Coupling Zone (LCZ). This LCZ can be related to tectonic structures, such as the Challenger Fracture Zone (CFZ) (Contreras-Reyes *et al.* 2015; Maksymowicz 2015). Poli *et al.* (2017) investigated the different fracture zones in the Illapel region (the CFZ, and the Juan Fernández Ridge, along with secondary structures), and suggested that these structures prevented the rupture to propagate further north and south. Consistently, we observe that the northern end of our coseismic slip zone correlates well with the CFZ. However, we do not find any correlation between
<figure><img src="image_23.png" /><figcaption>**Figure 14.** (a) Posterior mean coseismic slip model for a static inversion with a non shallow slip *a priori*. Arrows represent the slip directions and the ellipses their associated uncertainties. (b) Comparisons between tsunami observations (black) and stochastic predictions with shallow slip (red) and without shallow slip (blue). The tsunami waveform signal used in the inversion is shown between yellow dots. The map shows the depicted tsunami stations in (a).</figcaption></figure>

<table><caption>**Table 2.** BIC and AIC values with and without shallow slip. Bayesian (BIC) and Akaike (AIC) information criteria are defined in supplementary text S5. ΔBIC and ΔAIC are the difference in BIC and AIC values with respect to the slip model including shallow slip. The values suggest that the shallow slip should be included to properly explain the observations.</caption><thead><tr><td>Model</td><td>ΔBIC</td><td>ΔAIC</td></tr></thead><tbody><tr><td>Shallow slip (348 parameters)</td><td>0</td><td>0</td></tr><tr><td>No shallow slip (314 parameters)</td><td>1001</td><td>1096</td></tr></tbody></table>

since the maximum rise-time is around 30 s, which is considerably shorter than the total rupture time (around 100 s, Heaton 1990).

The encircled slip pulses visible in our solution between 30 and 60 s are consistent with previous back-projection results that suggest such complexities in the rupture (e.g. Ruiz *et al.* 2016; Meng *et al.* 2018). Ruiz *et al.* (2016) show an early stage bilateral rupture that later merged and propagated updip. Meng *et al.* (2018) report two episodes of splitting of rupture fronts, occurring both before reaching 60 s [an effect known as ‘double encircling pincer movement’ (Das & Kostrov 1983)]. The first episode reported by Meng *et al.* (2018) is between 15 and 35 s, and the second, around 45 and 60 s. The first encircling is colocated with the static coseismic model of An & Meng (2017). Consequently, Meng *et al.* (2018) suggest that the encircled region is an asperity. However, this static coseismic model could miss rupture features retrieved by our joint inversion that incorporates additional static and kinematic data. As previously pointed out (Ishii *et al.* 2007; Tilmann *et al.* 2016), back-projection sources trace the progression and changes of the rupture but are not proportional to slip. Our solution is more heterogenous, presenting multiple slip areas with both encircling episodes contouring

regions with small slip rates (and moderate slip), generating particularly high slip rates where the rupture focuses in the final stage of the earthquake (see time=60s, in Figure 16). In this sense, our observations suggest rather the contouring of two regions that do not slip during the rupture. Such strong changes in the rupture propagation associated with high slip rates explain the back-projection results of Meng *et al.* (2018). The small slip amplitude inside the contoured regions can be caused by different factors: (i) these areas could correspond to coupled regions (preventing seismic slip to occur), (ii) complexities at the subduction interface (e.g. due to fracture zones or seamounts) could prevent slip to propagate in these areas or (iii) the contoured regions could be far from the rupture (i.e. initial and dynamic stresses smaller than the fault strength). Regarding the coupling at the subduction interface, the model of the region proposed by Vigny *et al.* (2009) and updated by Métois *et al.* (2012) and Métois *et al.* (2016) shows a relatively high coupling coefficient in the Illapel earthquake area, except in the shallowest region, where the coefficient can be as low as 0.2. However, coupling close to the trench is usually poorly constrained by land-based geodetic data. The Illapel earthquake occurred in the Metropolitan segment defined by Métois *et al.* (2016), and is bounded in the north by the La Serena Low-Coupling Zone (LCZ). This LCZ can be related to tectonic structures, such as the Challenger Fracture Zone (CFZ) (Contreras-Reyes *et al.* 2015; Maksymowicz 2015). Poli *et al.* (2017) investigated the different fracture zones in the Illapel region (the CFZ, and the Juan Fernández Ridge, along with secondary structures), and suggested that these structures prevented the rupture to propagate further north and south. Consistently, we observe that the northern end of our coseismic slip zone correlates well with the CFZ. However, we do not find any correlation between
<figure><img src="image_24.png" /><figcaption>**Figure 15.** Comparison of posterior coseismic mean model with ISC aftershocks locations (green dots) after 12 hr (a), 24 hr (b) and 1 week after the main shock (c). The red star is the inverted hypocentre location. Arrows represent the slip directions with their corresponding uncertainty shown as ellipses.</figcaption></figure>

<figure><img src="image_25.png" /><figcaption>**Figure 16.** Five seconds snapshots of slip rate evolution. Slip rate is calculated using the posterior mean coseismic model considering the 2nd order $C_p$ solution. The red star is the inverted hypocentre location. Arrow lines represent the possible encircling locations.</figcaption></figure>

fault zone structures reported by Poli *et al.* (2017) and the encircled areas in our model. The small slip amplitudes in the contoured regions are thus likely not caused by such structures in the subducting plate.

To further investigate these encircled areas, we compare their locations with aftershocks distribution shown in Fig. 15. During the first 12 hr, we do not observe any aftershocks overlying the encircled regions. For the southern region, aftershocks depict a half
<figure><img src="image_24.png" /><figcaption>**Figure 15.** Comparison of posterior coseismic mean model with ISC aftershocks locations (green dots) after 12 hr (a), 24 hr (b) and 1 week after the main shock (c). The red star is the inverted hypocentre location. Arrows represent the slip directions with their corresponding uncertainty shown as ellipses.</figcaption></figure>

<figure><img src="image_25.png" /><figcaption>**Figure 16.** Five seconds snapshots of slip rate evolution. Slip rate is calculated using the posterior mean coseismic model considering the 2nd order $C_p$ solution. The red star is the inverted hypocentre location. Arrow lines represent the possible encircling locations.</figcaption></figure>

fault zone structures reported by Poli *et al.* (2017) and the encircled areas in our model. The small slip amplitudes in the contoured regions are thus likely not caused by such structures in the subducting plate.

To further investigate these encircled areas, we compare their locations with aftershocks distribution shown in Fig. 15. During the first 12 hr, we do not observe any aftershocks overlying the encircled regions. For the southern region, aftershocks depict a half
<figure><img src="image_26.png" /><figcaption>**Figure 17.** Posterior mean peak-slip rates. Slip rate is calculated using the posterior mean coseismic model using the 2nd order $C_p$ solution. Arrows represent the slip directions with their corresponding uncertainty. The red star is the inverted hypocentre location. Black contours show the posterior mean final slip model.</figcaption></figure>

semi-circle pattern that correlates well with our results. One week after the main shock, we note that both encircled regions remain with no significant aftershock activity. This is also shown in the two cross sections of Fig. A21, showing the absence of aftershocks in the encircled regions (shown as circles in Fig. A21). Several studies have linked aftershock occurrence with afterslip expansion over time (Kato 2007; Lengliné et al. 2012; Perfettini et al. 2018), often surrounding moderate/large coseismic slip areas (Mendoza & Hartzell 1988). Some fault areas around the Illapel rupture follow this behaviour, with an increase in the aftershocks rate, probably accompanying post-seismic slip in regions surrounding high co-seismic slip (cf. downdip slip region in Fig. 15a). However, the encircled areas remain seismically inactive after the main shock. The absence of aftershocks thus suggests that afterslip does not penetrate through these regions. Furthermore, according to the results of Frank et al. (2017), these two regions do not present any significant activity 9 months before, and 1 yr after the main shock. This suggests that the region would constitute a high-strength zone (i.e. with a high yield stress) compared with its surroundings (which could potentially break in the future), a region with a low slip deficit that broke recently (i.e. low initial stress), or with a larger fracture energy (Gallovič et al. 2020). The presence of high strength barriers has been observed for other megathrust earthquakes such as the 2001 $M_w$ = 8.1 Peru earthquake (Robinson et al. 2006), which was also associated with a low aftershock seismicity rate in the barrier region. On the other hand, if we consider the 1943 earthquake that occurred in the same region, and consider a fully coupled fault with a convergence rate of 67 mm yr$^{-1}$, the slip deficit would be 4.9 m, which is small compared to adjacent areas that experienced slip up to 20 m (cf. Fig. 6). If we take this slip deficit and calculate the corresponding scalar moment, we obtain a $M_0$ = 4.98 × 10$^{19}$ N·m ($M_w$ = 7.06) if they break individually, and $M_0$ = 4.48 × 10$^{20}$ N·m ($M_w$ = 7.7) if they break together.

## 6 CONCLUSION

Using extensive geodetic, seismic and tsunami data sets and a realistic uncertainty model, we obtain fully Bayesian finite-fault solutions of the 2015 $M_w$ = 8.3 Illapel earthquake. We employ a fixed subfault geometry and a non-linear parametrization (inverting for slip, rupture velocity, rise time and hypocentre location), which allows us to resolve the complexity of the rupture. We also propose a 2nd order perturbation approach to better account for prediction uncertainty in seismic waveforms.

Our kinematic slip models indicate two main slip asperities: a first asperity close to the hypocentre and another one at a shallow depth. Our analysis shows that shallow slip is required to fit tsunami observations and is consistent with the distribution of outer-rise aftershock seismicity. Historical records suggest that such shallow slip did not occur during the 1943 earthquake that affected the same region of the Chilean megathrust.

Our results also highlight encircling behaviours that occur when the rupture propagates towards the trench. Such rupture complexities have been previously suggested by back-projection studies. We suggest that these encircled regions are linked to areas associated with initial and dynamic stresses smaller than the fault yield stress. Further investigations are necessary to understand whether these areas correspond to low slip deficit regions or to fault areas with high strength that could be hosting future large earthquakes.

## ACKNOWLEDGMENTS

We thank H. Aochi and H. Bhat for helpful discussion. This project has received funding from the European Research Council (ERC, under the European Union's Horizon 2020 research and innovation programme under grant agreement no. 805256 and grant agreement no. 758210) and from Agence Nationale de la Recherche (project ANR-17-ERC3-0010). This research was also supported by the Mexican National Council for Science and Technology (CONACYT), scholarship 2018-000003-01EXTF-00012. RJ acknowledges funding from the Institut Universitaire de France. CL was supported by National Natural Science Foundation of China under grant 42274026. A portion of this work was conducted by the Jet Propulsion Laboratory, California Institute of Technology, under contract with the National Aeronautics and Space Administration. This work contains modified Copernicus data from the Sentinel-1A satellite processed by the ESA. We thank an anonymous reviewer, Frantisek Gallovič and Editor Eiichi Fukuyama for their valuable comments which improved this paper.

## DATA AVAILABILITY

The seismological data used in this study were acquired by CSN (Universidad de Chile 2012) and is freely accessible at the URL http://evtdb.csn.uchile.cl/. GPS displacements are available in Klein et al. (2017). Tide gauges are locally operated by Servicio Hidrográfico y Oceanográfico de la Armada and can be accessed at the URL http://www.ioc-sealevelmonitoring.org. The National Oceanic and Atmospheric Administration (NOAA) manages the DART stations accessible at https://www.ngdc.noaa.gov/hazard/DARTData.shtml (National Administration Oceanic and Atmospheric, 2005). InSAR images were acquired by the Sentinel-1A satellite operated by the European Space Agency under the Copernicus program and raw data can be consulted at https://winsar.unavco.org/data/access.
<figure><img src="image_26.png" /><figcaption>**Figure 17.** Posterior mean peak-slip rates. Slip rate is calculated using the posterior mean coseismic model using the 2nd order $C_p$ solution. Arrows represent the slip directions with their corresponding uncertainty. The red star is the inverted hypocentre location. Black contours show the posterior mean final slip model.</figcaption></figure>

semi-circle pattern that correlates well with our results. One week after the main shock, we note that both encircled regions remain with no significant aftershock activity. This is also shown in the two cross sections of Fig. A21, showing the absence of aftershocks in the encircled regions (shown as circles in Fig. A21). Several studies have linked aftershock occurrence with afterslip expansion over time (Kato 2007; Lengliné et al. 2012; Perfettini et al. 2018), often surrounding moderate/large coseismic slip areas (Mendoza & Hartzell 1988). Some fault areas around the Illapel rupture follow this behaviour, with an increase in the aftershocks rate, probably accompanying post-seismic slip in regions surrounding high co-seismic slip (cf. downdip slip region in Fig. 15a). However, the encircled areas remain seismically inactive after the main shock. The absence of aftershocks thus suggests that afterslip does not penetrate through these regions. Furthermore, according to the results of Frank et al. (2017), these two regions do not present any significant activity 9 months before, and 1 yr after the main shock. This suggests that the region would constitute a high-strength zone (i.e. with a high yield stress) compared with its surroundings (which could potentially break in the future), a region with a low slip deficit that broke recently (i.e. low initial stress), or with a larger fracture energy (Gallovič et al. 2020). The presence of high strength barriers has been observed for other megathrust earthquakes such as the 2001 $M_w$ = 8.1 Peru earthquake (Robinson et al. 2006), which was also associated with a low aftershock seismicity rate in the barrier region. On the other hand, if we consider the 1943 earthquake that occurred in the same region, and consider a fully coupled fault with a convergence rate of 67 mm yr$^{-1}$, the slip deficit would be 4.9 m, which is small compared to adjacent areas that experienced slip up to 20 m (cf. Fig. 6). If we take this slip deficit and calculate the corresponding scalar moment, we obtain a $M_0$ = 4.98 × 10$^{19}$ N·m ($M_w$ = 7.06) if they break individually, and $M_0$ = 4.48 × 10$^{20}$ N·m ($M_w$ = 7.7) if they break together.

## 6 CONCLUSION

Using extensive geodetic, seismic and tsunami data sets and a realistic uncertainty model, we obtain fully Bayesian finite-fault solutions of the 2015 $M_w$ = 8.3 Illapel earthquake. We employ a fixed subfault geometry and a non-linear parametrization (inverting for slip, rupture velocity, rise time and hypocentre location), which allows us to resolve the complexity of the rupture. We also propose a 2nd order perturbation approach to better account for prediction uncertainty in seismic waveforms.

Our kinematic slip models indicate two main slip asperities: a first asperity close to the hypocentre and another one at a shallow depth. Our analysis shows that shallow slip is required to fit tsunami observations and is consistent with the distribution of outer-rise aftershock seismicity. Historical records suggest that such shallow slip did not occur during the 1943 earthquake that affected the same region of the Chilean megathrust.

Our results also highlight encircling behaviours that occur when the rupture propagates towards the trench. Such rupture complexities have been previously suggested by back-projection studies. We suggest that these encircled regions are linked to areas associated with initial and dynamic stresses smaller than the fault yield stress. Further investigations are necessary to understand whether these areas correspond to low slip deficit regions or to fault areas with high strength that could be hosting future large earthquakes.

## ACKNOWLEDGMENTS

We thank H. Aochi and H. Bhat for helpful discussion. This project has received funding from the European Research Council (ERC, under the European Union's Horizon 2020 research and innovation programme under grant agreement no. 805256 and grant agreement no. 758210) and from Agence Nationale de la Recherche (project ANR-17-ERC3-0010). This research was also supported by the Mexican National Council for Science and Technology (CONACYT), scholarship 2018-000003-01EXTF-00012. RJ acknowledges funding from the Institut Universitaire de France. CL was supported by National Natural Science Foundation of China under grant 42274026. A portion of this work was conducted by the Jet Propulsion Laboratory, California Institute of Technology, under contract with the National Aeronautics and Space Administration. This work contains modified Copernicus data from the Sentinel-1A satellite processed by the ESA. We thank an anonymous reviewer, Frantisek Gallovič and Editor Eiichi Fukuyama for their valuable comments which improved this paper.

## DATA AVAILABILITY

The seismological data used in this study were acquired by CSN (Universidad de Chile 2012) and is freely accessible at the URL http://evtdb.csn.uchile.cl/. GPS displacements are available in Klein et al. (2017). Tide gauges are locally operated by Servicio Hidrográfica y Oceanográfico de la Armada and can be accessed at the URL http://www.ioc-sealevelmonitoring.org. The National Oceanic and Atmospheric Administration (NOAA) manages the DART stations accessible at https://www.ngdc.noaa.gov/hazard/DARTData.shtml (National Administration Oceanic and Atmospheric, 2005). InSAR images were acquired by the Sentinel-1A satellite operated by the European Space Agency under the Copernicus program and raw data can be consulted at https://winsar.unavco.org/data/access.
## SUPPORTING INFORMATION

Supplementary data are available at GJI online.

suppl.data

Please note: Oxford University Press is not responsible for the content or functionality of any supporting materials supplied by the authors. Any queries (other than missing material) should be directed to the corresponding author for the paper.

## REFERENCES

* National Oceanic and Atmospheric Administration. 2005. Deep-Ocean Assessment and Reporting of Tsunamis (DART(R)). NOAA National Centers for Environmental Information.

* An, C. & Meng, L., 2017. Time reversal imaging of the 2015 illapel tsunami source, *Geophys. Res. Lett.*, **44**(4), 1732-1739.

* An, C., Yue, H., Sun, J., Meng, L. & Báez, J.C., 2017. The 2015 mw 8.3 illapel, chile, earthquake: Direction-reversed along-dip rupture with localized water reverberation the 2015 Mw 8.3 Illapel, Chile, earthquake, *Bull. seism. Soc. Am.*, **107**(5), 2416-2426.

* Angermann, D., Klotz, J. & Reigber, C., 1999. Space-geodetic estimation of the Nazca-South America euler vector, *Earth planet. Sci. Lett.*, **171**(3), 329-334.

* Beck, S., Barrientos, S., Kausel, E. & Reyes, M., 1998. Source characteristics of historic earthquakes along the central Chile subduction Askew et Alzone, *J. South Am. Earth Sci.*, **11**(2), 115-129.

* Beresnev, I.A., 2003. Uncertainties in finite-fault slip inversions: to what extent to believe? (a critical review), *Bull. seism. Soc. Am.*, **93**(6), 2445-2458.

* Bishop, C.M., 2006. *Pattern Recognition and Machine Learning*, Springer.

* Bletery, Q., Sladen, A., Jiang, J. & Simons, M., 2016. A Bayesian source model for the 2004 great Sumatra-Andaman earthquake, *J. geophys. Res.*, **121**(7), 5116-5135.

* Cohee, B.P. & Beroza, G.C., 1994. A comparison of two methods for earthquake source inversion using strong motion seismograms, *Ann. Geophys.*, XXXVII(6), doi:10.4401/ag-4151.

* Contreras-Reyes, E., Ruiz, J.A., Becerra, J., Kopp, H., Reichert, C., Maksymowicz, A. & Arriagada, C., 2015. Structure and tectonics of the central chilean margin (31°-33° s): implications for subduction erosion and shallow crustal seismicity, *Geophys. J. Int.*, **203**(2), 776-791.

* Das, S. & Kostrov, B., 1983. Breaking of a single asperity: rupture process and seismic radiation, *J. geophys. Res.*, **88**(B5), 4277-4288.

* Duputel, Z. et al., 2015. The iquique earthquake sequence of April 2014: Bayesian modeling accounting for prediction uncertainty, *Geophys. Res. Lett.*, **42**(19), 7949-7957.

* Duputel, Z., Agram, P.S., Simons, M., Minson, S.E. & Beck, J.L., 2014. Accounting for prediction uncertainty when inferring subsurface fault slip, *Geophys. J. Int.*, **197**(1), 464-482.

* Duputel, Z., Rivera, L., Fukahata, Y. & Kanamori, H., 2012. Uncertainty estimations for seismic source inversions, *Geophys. J. Int.*, **190**(2), 1243-1256.

* 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, *J. geophys. Res.*, **86**(B4), 2825-2852.

* Ekström, G., Nettles, M. & Dziewoński, A., 2012. The global CMT project 2004-2010: centroid-moment tensors for 13,017 earthquakes, *Phys. Earth planet. Inter.*, **200**, 1-9.

* Fernández, J., Pastén, C., Ruiz, S. & Leyton, F., 2019. Damage assessment of the 2015 Mw 8.3 Illapel earthquake in the North-Central Chile, *Nat. Hazards*, **96**(1), 269-283.

* Frank, W.B., Poli, P. & Perfettini, H., 2017. Mapping the rheology of the central Chile subduction zone with aftershocks, *Geophys. Res. Lett.*, **44**(11), 5374-5382.

* Gallovič, F., Zahradník, J., Plicka, V., Sokos, E., Evangelidis, C., Fountoulakis, I. & Turhan, F., 2020. Complex rupture dynamics on an immature fault during the 2020 Mw 6.8 Elazığ earthquake, Turkey, *Commun. Earth Environ.*, **1**(1), 40.

* Gombert, B., Duputel, Z., Jolivet, R., Doubre, C., Rivera, L. & Simons, M., 2018a. Revisiting the 1992 landers earthquake: a Bayesian exploration of co-seismic slip and off-fault damage, *Geophys. J. Int.*, **212**(2), 839-852.

* Gombert, B., Duputel, Z., Jolivet, R., Simons, M., Jiang, J., Liang, C., Fielding, E.J. & Rivera, L., 2018b. Strain budget of the Ecuador-Colombia subduction zone: a stochastic view, *Earth planet. Sci. Lett.*, **498**, 288-299.

* Hallo, M. & Gallovič, F., 2016. Fast and cheap approximation of Green's function uncertainty for waveform-based earthquake source inversions, *Geophys. J. Int.*, **207**(2), 1012-1029.

* Hartzell, S.H. & Heaton, T.H., 1983. Inversion of strong ground motion and teleseismic waveform data for the fault rupture history of the 1979 Imperial Valley, California, earthquake, *Bull. seism. Soc. Am.*, **73**(6A), 1553-1583.

* Heaton, T.H., 1990. Evidence for and implications of self-healing pulses of slip in earthquake rupture, *Phys. Earth planet. Inter.*, **64**(1), 1-20.

* Heidarzadeh, M., Murotani, S., Satake, K., Ishibe, T. & Gusman, A.R., 2016. Source model of the 16 September 2015 Illapel, Chile, Mw 8.4 earthquake based on teleseismic and tsunami data, *Geophys. Res. Lett.*, **43**(2), 643-650.

* Herrmann, R.B., 2013. Computer programs in seismology: an evolving tool for instruction and research, *Seismol. Res. Lett.*, **84**(6), 1081-1088.

* Ide, S., 2015. 4.09-Slip inversion, *Treatise on Geophysics*, 2nd edn, pp. 215-241, Elsevier.

* Ishii, M., Shearer, P.M., Houston, H. & Vidale, J.E., 2007. Teleseismic P wave imaging of the 26 December 2004 Sumatra-Andaman and 28 March 2005 Sumatra earthquake ruptures using the hi-net array, *J. geophys. Res.*, **112**(B11), doi:10.1029/2006JB004700.

* Jolivet, R., Jara, J., Dalaison, M., Rouet-Leduc, B., Özdemir, A., Dogan, U., Çakir, Z. & Ergintav, S., 2023. Daily to centennial behavior of aseismic slip along the central section of the North Anatolian Fault, *J. geophys. Res.*, **128**(7), e2022JB026018, doi:10.1029/2022JB026018.

* Jolivet, R., Simons, M., Agram, P., Duputel, Z. & Shen, Z.-K., 2015. Aseismic slip and seismogenic coupling along the Central San Andreas Fault, *Geophys. Res. Lett.*, **42**(2), 297-306.

* Jolivet, R., Simons, M., Duputel, Z., Olive, J.-A., Bhat, H. & Bletery, Q., 2020. Interseismic loading of subduction megathrust drives long-term uplift in Northern Chile, *Geophys. Res. Lett.*, **47**(8), e2019GL085377, doi:10.1029/2019GL085377.

* Kato, N., 2007. Expansion of aftershock areas caused by propagating post-seismic sliding, *Geophys. J. Int.*, **168**(2), 797-808.

* 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 planet. Sci. Lett.*, **469**, 123-134.

* 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, Mw 8.3 earthquake, *Geophys. Res. Lett.*, **43**(15), 7902-7912.

* Lee, S.-J., Yeh, T.-Y., Lin, T.-C., Lin, Y.-Y., Song, T.-R.A. & Huang, B.-S., 2016. Two-stage composite megathrust rupture of the 2015 Mw8.4 Illapel, Chile, earthquake identified by spectral-element inversion of teleseismic waves, *Geophys. Res. Lett.*, **43**(10), 4979-4985.

* Lengliné, O., Enescu, B., Peng, Z. & Shiomi, K., 2012. Decay and expansion of the early aftershock activity following the 2011, Mw9.0 Tohoku earthquake, *Geophys. Res. Lett.*, **39**(18), doi:10.1029/2012GL052797.

* 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, Mw 8.3 earthquake rupture process, *Geophys. Res. Lett.*, **43**(9), 4303-4312.

* Lindsey, E.O., Mallick, R., Hubbard, J.A., Bradley, K.E., Almeida, R.V., Moore, J.D., Bürgmann, R. & Hill, E.M., 2021. Slip rate deficit and earthquake potential on shallow megathrusts, *Nat. Geosci.*, **14**(5), 321-326.
## SUPPORTING INFORMATION

Supplementary data are available at GJI online.

### suppl.data

Please note: Oxford University Press is not responsible for the content or functionality of any supporting materials supplied by the authors. Any queries (other than missing material) should be directed to the corresponding author for the paper.

## REFERENCES

National Oceanic and Atmospheric Administration. 2005. Deep-Ocean Assessment and Reporting of Tsunamis (DART(R)). NOAA National Centers for Environmental Information.

An, C. & Meng, L., 2017. Time reversal imaging of the 2015 illapel tsunami source, *Geophys. Res. Lett.*, **44**(4), 1732–1739.

An, C., Yue, H., Sun, J., Meng, L. & Báez, J.C., 2017. The 2015 mw 8.3 illapel, chile, earthquake: Direction-reversed along-dip rupture with localized water reverberation the 2015 Mw 8.3 Illapel, Chile, earthquake, *Bull. seism. Soc. Am.*, **107**(5), 2416–2426.

Angermann, D., Klotz, J. & Reigber, C., 1999. Space-geodetic estimation of the Nazca-South America euler vector, *Earth planet. Sci. Lett.*, **171**(3), 329–334.

Beck, S., Barrientos, S., Kausel, E. & Reyes, M., 1998. Source characteristics of historic earthquakes along the central Chile subduction Askew et Alzone, *J. South Am. Earth Sci.*, **11**(2), 115–129.

Beresnev, I.A., 2003. Uncertainties in finite-fault slip inversions: to what extent to believe? (a critical review), *Bull. seism. Soc. Am.*, **93**(6), 2445–2458.

Bishop, C.M., 2006. *Pattern Recognition and Machine Learning*, Springer.

Bletery, Q., Sladen, A., Jiang, J. & Simons, M., 2016. A Bayesian source model for the 2004 great Sumatra-Andaman earthquake, *J. geophys. Res.*, **121**(7), 5116–5135.

Cohee, B.P. & Beroza, G.C., 1994. A comparison of two methods for earthquake source inversion using strong motion seismograms, *Ann. Geophys.*, **XXXVII**(6), doi:10.4401/ag-4151.

Contreras-Reyes, E., Ruiz, J.A., Becerra, J., Kopp, H., Reichert, C., Maksymowicz, A. & Arriagada, C., 2015. Structure and tectonics of the central chilean margin (31°–33° s): implications for subduction erosion and shallow crustal seismicity, *Geophys. J. Int.*, **203**(2), 776–791.

Das, S. & Kostrov, B., 1983. Breaking of a single asperity: rupture process and seismic radiation, *J. geophys. Res.*, **88**(B5), 4277–4288.

Duputel, Z. et al., 2015. The iquique earthquake sequence of April 2014: Bayesian modeling accounting for prediction uncertainty, *Geophys. Res. Lett.*, **42**(19), 7949–7957.

Duputel, Z., Agram, P.S., Simons, M., Minson, S.E. & Beck, J.L., 2014. Accounting for prediction uncertainty when inferring subsurface fault slip, *Geophys. J. Int.*, **197**(1), 464–482.

Duputel, Z., Rivera, L., Fukahata, Y. & Kanamori, H., 2012. Uncertainty estimations for seismic source inversions, *Geophys. J. Int.*, **190**(2), 1243–1256.

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, *J. geophys. Res.*, **86**(B4), 2825–2852.

Ekström, G., Nettles, M. & Dziewoński, A., 2012. The global CMT project 2004–2010: centroid-moment tensors for 13,017 earthquakes, *Phys. Earth planet. Inter.*, **200**, 1–9.

Fernández, J., Pastén, C., Ruiz, S. & Leyton, F., 2019. Damage assessment of the 2015 Mw 8.3 Illapel earthquake in the North-Central Chile, *Nat. Hazards*, **96**(1), 269–283.

Frank, W.B., Poli, P. & Perfettini, H., 2017. Mapping the rheology of the central Chile subduction zone with aftershocks, *Geophys. Res. Lett.*, **44**(11), 5374–5382.

Gallovič, F., Zahradník, J., Plicka, V., Sokos, E., Evangelidis, C., Fountoulakis, I. & Turhan, F., 2020. Complex rupture dynamics on an immature fault during the 2020 Mw 6.8 Elazığ earthquake, Turkey, *Commun. Earth Environ.*, **1**(1), 40.

Gombert, B., Duputel, Z., Jolivet, R., Doubre, C., Rivera, L. & Simons, M., 2018a. Revisiting the 1992 landers earthquake: a Bayesian exploration of co-seismic slip and off-fault damage, *Geophys. J. Int.*, **212**(2), 839–852.

Gombert, B., Duputel, Z., Jolivet, R., Simons, M., Jiang, J., Liang, C., Fielding, E.J. & Rivera, L., 2018b. Strain budget of the Ecuador-Colombia subduction zone: a stochastic view, *Earth planet. Sci. Lett.*, **498**, 288–299.

Hallo, M. & Gallovič, F., 2016. Fast and cheap approximation of Green's function uncertainty for waveform-based earthquake source inversions, *Geophys. J. Int.*, **207**(2), 1012–1029.

Hartzell, S.H. & Heaton, T.H., 1983. Inversion of strong ground motion and teleseismic waveform data for the fault rupture history of the 1979 Imperial Valley, California, earthquake, *Bull. seism. Soc. Am.*, **73**(6A), 1553–1583.

Heaton, T.H., 1990. Evidence for and implications of self-healing pulses of slip in earthquake rupture, *Phys. Earth planet. Inter.*, **64**(1), 1–20.

Heidarzadeh, M., Murotani, S., Satake, K., Ishibe, T. & Gusman, A.R., 2016. Source model of the 16 September 2015 Illapel, Chile, Mw 8.4 earthquake based on teleseismic and tsunami data, *Geophys. Res. Lett.*, **43**(2), 643–650.

Herrmann, R.B., 2013. Computer programs in seismology: an evolving tool for instruction and research, *Seismol. Res. Lett.*, **84**(6), 1081–1088.

Ide, S., 2015. 4.09-Slip inversion, *Treatise on Geophysics*, 2nd edn, pp. 215–241, Elsevier.

Ishii, M., Shearer, P.M., Houston, H. & Vidale, J.E., 2007. Teleseismic P wave imaging of the 26 December 2004 Sumatra-Andaman and 28 March 2005 Sumatra earthquake ruptures using the hi-net array, *J. geophys. Res.*, **112**(B11), doi:10.1029/2006JB004700.

Jolivet, R., Jara, J., Dalaison, M., Rouet-Leduc, B., Özdemir, A., Dogan, U., Çakir, Z. & Ergintav, S., 2023. Daily to centennial behavior of aseismic slip along the central section of the North Anatolian Fault, *J. geophys. Res.*, **128**(7), e2022JB026018, doi:10.1029/2022JB026018.

Jolivet, R., Simons, M., Agram, P., Duputel, Z. & Shen, Z.-K., 2015. Aseismic slip and seismogenic coupling along the Central San Andreas Fault, *Geophys. Res. Lett.*, **42**(2), 297–306.

Jolivet, R., Simons, M., Duputel, Z., Olive, J.-A., Bhat, H. & Bletery, Q., 2020. Interseismic loading of subduction megathrust drives long-term uplift in Northern Chile, *Geophys. Res. Lett.*, **47**(8), e2019GL085377, doi:10.1029/2019GL085377.

Kato, N., 2007. Expansion of aftershock areas caused by propagating post-seismic sliding, *Geophys. J. Int.*, **168**(2), 797–808.

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 planet. Sci. Lett.*, **469**, 123–134.

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, Mw 8.3 earthquake, *Geophys. Res. Lett.*, **43**(15), 7902–7912.

Lee, S.-J., Yeh, T.-Y., Lin, T.-C., Lin, Y.-Y., Song, T.-R.A. & Huang, B.-S., 2016. Two-stage composite megathrust rupture of the 2015 Mw8.4 Illapel, Chile, earthquake identified by spectral-element inversion of teleseismic waves, *Geophys. Res. Lett.*, **43**(10), 4979–4985.

Lengliné, O., Enescu, B., Peng, Z. & Shiomi, K., 2012. Decay and expansion of the early aftershock activity following the 2011, Mw9.0 Tohoku earthquake, *Geophys. Res. Lett.*, **39**(18), doi:10.1029/2012GL052797.

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, Mw 8.3 earthquake rupture process, *Geophys. Res. Lett.*, **43**(9), 4303–4312.

Lindsey, E.O., Mallick, R., Hubbard, J.A., Bradley, K.E., Almeida, R.V., Moore, J.D., Bürgmann, R. & Hill, E.M., 2021. Slip rate deficit and earthquake potential on shallow megathrusts, *Nat. Geosci.*, **14**(5), 321–326.
Lomnitz, C., 2004. Major earthquakes of Chile: a historical survey, 1535-1960, *Seismol. Res. Lett.*, **75**(3), 368-378.

Maksymowicz, A., 2015. The geometry of the Chilean continental wedge: tectonic segmentation of subduction processes off chile, *Tectonophysics*, **659**, 183-196.

Melgar, D. et al., 2016. Slip segmentation and slow rupture to the trench during the 2015, Mw8. 3 Illapel, Chile earthquake, *Geophys. Res. Lett.*, **43**(3), 961-966.

Mendoza, C. & Hartzell, S.H., 1988. Aftershock patterns and main shock faulting, *Bull. seism. Soc. Am.*, **78**(4), 1438-1449.

Meng, L., Bao, H., Huang, H., Zhang, A., Bloore, A. & Liu, Z., 2018. Double pincer movement: encircling rupture splitting during the 2015 Mw 8.3 Illapel earthquake, *Earth planet. Sci. Lett.*, **495**, 164-173.

Métois, M., Socquet, A. & Vigny, C., 2012. Interseismic coupling, segmentation and mechanical behavior of the central Chile subduction zone, *J. geophys. Res.*, **117**(B3), doi:10.1029/2011JB008736.

Métois, M., Vigny, C. & Socquet, A., 2016. Interseismic coupling, megathrust earthquakes and seismic swarms along the Chilean subduction zone (38-18 s), *Pure appl. Geophys.*, **173**(5), 1431-1449.

Minson, S., Simons, M. & Beck, J., 2013. Bayesian inversion for finite fault earthquake source models I—Theory and algorithm, *Geophys. J. Int.*, **194**(3), 1701-1726.

Okuwaki, R., Yagi, Y., Aránguiz, R., González, J. & González, G., 2016. Rupture process during the 2015 Illapel, Chile earthquake: zigzag-along-dip rupture episodes, *Pure appl. Geophys.*, **173**(4), 1011-1020.

Page, M.T., Dunham, E.M. & Carlson, J., 2005. Distinguishing barriers and asperities in near-source ground motion, *J. geophys. Res.*, **110**(B11), doi:10.1029/2005JB003736.

Perfettini, H., Frank, W., Marsan, D. & Bouchon, M., 2018. A model of aftershock migration driven by afterslip, *Geophys. Res. Lett.*, **45**(5), 2283-2293.

Poli, P., Maksymowicz, A. & Ruiz, S., 2017. The Mw 8.3 Illapel earthquake (Chile): preseismic and postseismic activity associated with hydrated slab structures, *Geology*, **45**(3), 247-250.

Ragon, T., Sladen, A. & Simons, M., 2018. Accounting for uncertain fault geometry in earthquake source inversions–I. Theory and simplified application, *Geophys. J. Int.*, **214**(2), 1174-1190.

Razafindrakoto, H.N. & Mai, P.M., 2014. Uncertainty in earthquake source imaging due to variations in source time function and Earth structure, *Bull. seism. Soc. Am.*, **104**(2), 855-874.

Robinson, D., Das, S. & Watts, A., 2006. Earthquake rupture stalled by a subducting fracture zone, *Science*, **312**(5777), 1203-1205.

Ruiz, S. et al., 2016. The seismic sequence of the 16 September 2015 Mw 8.3 Illapel, Chile, earthquake, *Seismol. Res. Lett.*, **87**(4), 789-799.

Ruiz, S. & Madariaga, R., 2018. Historical and recent large megathrust earthquakes in Chile, *Tectonophysics*, **733**, 37-56.

Satake, K. & Heidarzadeh, M., 2017. A review of source models of the 2015 Illapel, Chile earthquake and insights from tsunami data, *Pure appl. Geophys.*, **174**(1), 1-9.

Shrivastava, M.N. et al., 2016. Coseismic slip and afterslip of the 2015 Mw 8.3 Illapel (Chile) earthquake determined from continuous GPS data, *Geophys. Res. Lett.*, **43**(20), 10-710.

Sladen, A. & Trevisan, J., 2018. Shallow megathrust earthquake ruptures betrayed by their outer-trench aftershocks signature, *Earth planet. Sci. Lett.*, **483**, 105-113.

Tilmann, F. et al., 2016. The 2015 Illapel earthquake, central Chile: a type case for a characteristic earthquake?, *Geophys. Res. Lett.*, **43**(2), 574-583.

Tromp, J., Tape, C. & Liu, Q., 2005. Seismic tomography, adjoint methods, time reversal and banana-doughnut kernels, *Geophys. J. Int.*, **160**(1), 195-216.

Twardzik, C., Vergnolle, M., Sladen, A. & Tsang, L.L., 2021. Very early identification of a bimodal frictional behavior during the post-seismic phase of the 2015 Mw 8.3 Illapel, Chile, earthquake, *Solid Earth*, **12**(11), 2523-2537.

Universidad de Chile, 2012. Red Sismologica Nacional [Data set]. International Federation of Digital Seismograph Networks. doi:10.7914/SN/C1.

Vigny, C., Rudloff, A., Ruegg, J.-C., Madariaga, R., Campos, J. & Alvarez, M., 2009. Upper plate deformation measured by GPS in the Coquimbo Gap, Chile, *Phys. Earth planet. Inter.*, **175**(1-2), 86-95.

Wald, D.J. & Graves, R.W., 2001. Resolution analysis of finite fault source inversion using one- and three-dimensional Green's functions: 2. Combining seismic and geodetic data, *J. geophys. Res.*, **106**(B5), 8767-8788.

Williams, C.A. & Wallace, L.M., 2015. Effects of material property variations on slip estimates for subduction interface slow-slip events, *Geophys. Res. Lett.*, **42**(4), 1113-1121.

Williamson, A., Newman, A. & Cummins, P., 2017. Reconstruction of co-seismic slip from the 2015 Illapel earthquake using combined geodetic and tsunami waveform data, *J. geophys. Res.*, **122**(3), 2119-2130.

Yagi, Y. & Fukahata, Y., 2011. Introduction of uncertainty of Green's function into waveform inversion for seismic source processes, *Geophys. J. Int.*, **186**(2), 711-720.

Yin, J., Yang, H., Yao, H. & Weng, H., 2016. Coseismic radiation and stress drop during the 2015 Mw 8.3 Illapel, Chile megathrust earthquake, *Geophys. Res. Lett.*, **43**(4), 1520-1528.

Zhu, L. & Rivera, L.A., 2002. A note on the dynamic and static displacements from a point source in multilayered media, *Geophys. J. Int.*, **148**(3), 619-627.
Lomnitz, C., 2004. Major earthquakes of Chile: a historical survey, 1535-1960, *Seismol. Res. Lett.*, **75**(3), 368-378.

Maksymowicz, A., 2015. The geometry of the Chilean continental wedge: tectonic segmentation of subduction processes off chile, *Tectonophysics*, **659**, 183-196.

Melgar, D. et al., 2016. Slip segmentation and slow rupture to the trench during the 2015, Mw8. 3 Illapel, Chile earthquake, *Geophys. Res. Lett.*, **43**(3), 961-966.

Mendoza, C. & Hartzell, S.H., 1988. Aftershock patterns and main shock faulting, *Bull. seism. Soc. Am.*, **78**(4), 1438-1449.

Meng, L., Bao, H., Huang, H., Zhang, A., Bloore, A. & Liu, Z., 2018. Double pincer movement: encircling rupture splitting during the 2015 Mw 8.3 Illapel earthquake, *Earth planet. Sci. Lett.*, **495**, 164-173.

Métois, M., Socquet, A. & Vigny, C., 2012. Interseismic coupling, segmentation and mechanical behavior of the central Chile subduction zone, *J. geophys. Res.*, **117**(B3), doi:10.1029/2011JB008736.

Métois, M., Vigny, C. & Socquet, A., 2016. Interseismic coupling, megathrust earthquakes and seismic swarms along the Chilean subduction zone (38-18 s), *Pure appl. Geophys.*, **173**(5), 1431-1449.

Minson, S., Simons, M. & Beck, J., 2013. Bayesian inversion for finite fault earthquake source models I—Theory and algorithm, *Geophys. J. Int.*, **194**(3), 1701-1726.

Okuwaki, R., Yagi, Y., Aránguiz, R., González, J. & González, G., 2016. Rupture process during the 2015 Illapel, Chile earthquake: zigzag-along-dip rupture episodes, *Pure appl. Geophys.*, **173**(4), 1011-1020.

Page, M.T., Dunham, E.M. & Carlson, J., 2005. Distinguishing barriers and asperities in near-source ground motion, *J. geophys. Res.*, **110**(B11), doi:10.1029/2005JB003736.

Perfettini, H., Frank, W., Marsan, D. & Bouchon, M., 2018. A model of aftershock migration driven by afterslip, *Geophys. Res. Lett.*, **45**(5), 2283-2293.

Poli, P., Maksymowicz, A. & Ruiz, S., 2017. The Mw 8.3 Illapel earthquake (Chile): preseismic and postseismic activity associated with hydrated slab structures, *Geology*, **45**(3), 247-250.

Ragon, T., Sladen, A. & Simons, M., 2018. Accounting for uncertain fault geometry in earthquake source inversions–I. Theory and simplified application, *Geophys. J. Int.*, **214**(2), 1174-1190.

Razafindrakoto, H.N. & Mai, P.M., 2014. Uncertainty in earthquake source imaging due to variations in source time function and Earth structure, *Bull. seism. Soc. Am.*, **104**(2), 855-874.

Robinson, D., Das, S. & Watts, A., 2006. Earthquake rupture stalled by a subducting fracture zone, *Science*, **312**(5777), 1203-1205.

Ruiz, S. et al., 2016. The seismic sequence of the 16 September 2015 Mw 8.3 Illapel, Chile, earthquake, *Seismol. Res. Lett.*, **87**(4), 789-799.

Ruiz, S. & Madariaga, R., 2018. Historical and recent large megathrust earthquakes in Chile, *Tectonophysics*, **733**, 37-56.

Satake, K. & Heidarzadeh, M., 2017. A review of source models of the 2015 Illapel, Chile earthquake and insights from tsunami data, *Pure appl. Geophys.*, **174**(1), 1-9.

Shrivastava, M.N. et al., 2016. Coseismic slip and afterslip of the 2015 Mw 8.3 Illapel (Chile) earthquake determined from continuous GPS data, *Geophys. Res. Lett.*, **43**(20), 10-710.

Sladen, A. & Trevisan, J., 2018. Shallow megathrust earthquake ruptures betrayed by their outer-trench aftershocks signature, *Earth planet. Sci. Lett.*, **483**, 105-113.

Tilmann, F. et al., 2016. The 2015 Illapel earthquake, central Chile: a type case for a characteristic earthquake?, *Geophys. Res. Lett.*, **43**(2), 574-583.

Tromp, J., Tape, C. & Liu, Q., 2005. Seismic tomography, adjoint methods, time reversal and banana-doughnut kernels, *Geophys. J. Int.*, **160**(1), 195-216.

Twardzik, C., Vergnolle, M., Sladen, A. & Tsang, L.L., 2021. Very early identification of a bimodal frictional behavior during the post-seismic phase of the 2015 Mw 8.3 Illapel, Chile, earthquake, *Solid Earth*, **12**(11), 2523-2537.

Universidad de Chile, 2012. Red Sismologica Nacional [Data set]. International Federation of Digital Seismograph Networks. doi:10.7914/SN/C1.

Vigny, C., Rudloff, A., Ruegg, J.-C., Madariaga, R., Campos, J. & Alvarez, M., 2009. Upper plate deformation measured by GPS in the Coquimbo Gap, Chile, *Phys. Earth planet. Inter.*, **175**(1-2), 86-95.

Wald, D.J. & Graves, R.W., 2001. Resolution analysis of finite fault source inversion using one- and three-dimensional Green's functions: 2. Combining seismic and geodetic data, *J. geophys. Res.*, **106**(B5), 8767-8788.

Williams, C.A. & Wallace, L.M., 2015. Effects of material property variations on slip estimates for subduction interface slow-slip events, *Geophys. Res. Lett.*, **42**(4), 1113-1121.

Williamson, A., Newman, A. & Cummins, P., 2017. Reconstruction of co-seismic slip from the 2015 Illapel earthquake using combined geodetic and tsunami waveform data, *J. geophys. Res.*, **122**(3), 2119-2130.

Yagi, Y. & Fukahata, Y., 2011. Introduction of uncertainty of Green's function into waveform inversion for seismic source processes, *Geophys. J. Int.*, **186**(2), 711-720.

Yin, J., Yang, H., Yao, H. & Weng, H., 2016. Coseismic radiation and stress drop during the 2015 Mw 8.3 Illapel, Chile megathrust earthquake, *Geophys. Res. Lett.*, **43**(4), 1520-1528.

Zhu, L. & Rivera, L.A., 2002. A note on the dynamic and static displacements from a point source in multilayered media, *Geophys. J. Int.*, **148**(3), 619-627.