# Fully joint inversion of the 2016 M w 7.6 Chiloé earthquake

F Bravo, S Peyrat, F Delgado, M Fuentes, B Derode, A Perez, J Campos

► To cite this version:

F Bravo, S Peyrat, F Delgado, M Fuentes, B Derode, et al.. Fully joint inversion of the 2016 M w 7.6 Chiloé earthquake. Geophysical Journal International, 2023, 232 (3), pp.2001-2016. 10.1093/gji/ggac411. hal-04727024

HAL Id: hal-04727024

https://hal.science/hal-04727024v1

Submitted on 9 Oct 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_1.png" /></figure>

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

# Fully joint inversion of the 2016 M<sub>w</sub> 7.6 Chiloé earthquake

F. Bravo<sup>1</sup>, S. Peyrat,<sup>2</sup> F. Delgado,<sup>3</sup> M. Fuentes,<sup>4</sup> B. Derode,<sup>4</sup> A. Perez<sup>4</sup> and J. Campos<sup>4</sup>

Accepted 2022 October 17. Received 2022 October 12; in original form 2021 September 8

## SUMMARY

A large (M<sub>w</sub> 7.6) megathrust earthquake occurred on 2016 December 25 in Southern Chile, south of the Chiloé Island (74.2°W, 43.3°S) in the South America–Nazca subduction zone. This earthquake was the first large event in this seismotectonic segment since the M<sub>w</sub> 9.5 1960 Valdivia megathrust earthquake and broke a ~50-km-long segment of the southern part of its rupture zone. Source parameters are inferred from teleseismic broad-bands, strong motions, GPS, cGPS, InSAR and tide gauge data. We show that the joint inversion significantly improves the resolution of the slip distribution, taking advantage of each data set. Our slip models predict a single slip patch of 70 km × 60 km with a maximum slip of 3.2 ± 0.8 m and a moment magnitude of 7.64. The hypothesis from previous studies that the Chiloé earthquake released energy accumulated before the Valdivia earthquake is not supported by the ensemble of geodetic, seismological and tsunami data. Hence, the Chiloé earthquake most likely released all the strain accumulated in the rupture area since the 1960 earthquake.

**Key words:** South America; Joint inversion; Earthquake source observations; Subduction zone processes.

## 1 INTRODUCTION

A substantial limitation in seismic risk assessment along subduction zones is due to the short windows of observations of a complete cycle of interseismic energy accumulation and release. Deducing the whole history of seismic cycle represents a critical goal to improve the potential seismic estimations along the high-populated continental margins. Correlations between the interseismic coupling and coseismic slip models show changes in plate locking during the entire seismic cycle. Precise determination of these changes opens the opportunity to improve the time-depended probability estimations of both the recurrence of the next great earthquake and its magnitude (Melnick et al. 2018).

During the last three decades, the megathrust earthquake toolkit has expanded from seismology to include space-geodetic and tide gauge data to better understand the distribution and correlation between coseismic and interseismic slips. Coseismic slip models are commonly estimated by inverting single or multiple data sources, including strong motion, teleseismic, GNSS, InSAR and tide gauge records. Further, retrospective studies have benefitted from the use of measurements that were not available at the time of the events (Yoshida & Koketsu 1990; Satake & Kanamori 1991; Wald & Heaton 1994). The main methodological problems incorporating data sets that measure a specific aspect of the earthquake rupture process are assessing the resolving power of each one with different spatial resolution restrictions, and the lack of offshore data close to

the trench. Therefore, a joint inversion of multiple data sets with different spatial and temporal resolution has the potential to offer more robust constrains on fault slip. However, both the availability of multiple data sets and adequate numerical tools that can handle them result in few examples of fully joint inversions (Yokota et al. 2011; Yue et al. 2014; Duputel et al. 2015).

On the other hand, plate coupling dynamics during the strain build-up in the interseismic-phase that follows a megathrust earthquake remains poorly understood. That is mainly due to the lack of decade-long geodetic records with enough spatial resolution of the plate coupling zone that record both the interseismic strain and the effect produced by the post-seismic viscoelastic relaxation that can last for several decades (Meade et al. 2013). One of the areas with a lack of detailed observations is the segment that broke during the Valdivia M<sub>w</sub> 9.5 megathrust in the Nazca-South America subduction zone. This earthquake is the largest subduction event recorded to date (e.g. Plafker & Savage 1970; Kanamori 1977; Cifuentes 1989; Barrientos & Ward 1990; Moreno et al. 2009; Ojeda et al. 2020). This giant earthquake activated the interplate contact along the southern part of the plate margin, along a trench-parallel segment of ~1000 km.

The inversion of terrestrial geodetic data of the 1960 megathrust showed a coseismic slip pattern characterized by four coseismic asperities (Fig. 1). The two largest asperities were located in front of the South American continent, while the two minor asperities were situated on the fjord-zone coast of southern Chile

<footnote n="1"><sup>1</sup>Centro Sismológico Nacional, Facultad de Ciencias Físicas y Matemáticas, Universidad de Chile, 8370415 Santiago, Chile. E-mail: fbravo@csn.uchile.cl</footnote>

<footnote n="2"><sup>2</sup>Géosciences Montpellier, Université de Montpellier, CNRS, 34095 Montpellier, France</footnote>

<footnote n="3"><sup>3</sup>Departamento de Geología, Centro de Excelencia en Geotermia de los Andes (CEGA), Universidad de Chile, 8370415 Santiago, Chile</footnote>

<footnote n="4"><sup>4</sup>Departamento de Geofísica, Facultad de Ciencias Físicas y Matemáticas, Universidad de Chile, 8370415 Santiago, Chile</footnote>
<figure><img src="image_3.png" /><figcaption>**Figure 1.** Tectonic setting: major structures are shown (Trench, LOFZ, Ridge). The fault traces showing the LOFZ are from the active and potentially active faults database of Maldonado *et al.* (2021). Seismicity from the International Seismological Centre Bulletin (1900 up to now). Yellow star is the epicentre of 2016 December 25 Chiloé earthquake. Blue star is the epicentre of the 1960 Valdivia earthquake. The isoslip curves (in m) show the rupture area of the 1960 Valdivia earthquake is from Moreno *et al.* (2009).</figcaption></figure>

(Moreno *et al.* 2011). The end of the rupture zone reaches the triple junction between Nazca–Antarctica and South America plates, which acts as a thermomechanical barrier that arrested the rupture (Contreras-Reyes & Carrizo 2011). The third relevant coseismic asperity is located immediately to the southwest of the Chiloé Island and east to Guafo Island. The southern asperity was located at the west of the Las Guaitecas Archipielago (Moreno *et al.* 2011). Despite the large magnitude of the earthquake, the seismological record shows low seismicity during the post-seismic phase that transition to interseismic strain build-up. The inversion of space geodetic data shows a heterogeneous interplate locking pattern (Moreno *et al.* 2011). In detail, an interseismic asperity that has a high-coupling is

located in the south of Chiloé Island, and fits well with the coseismic asperity located close to Guafo Island.

56 years after the 1960 $M_w$ 9.5 Valdivia earthquake, an offshore $M_w$ 7.6 earthquake occurred at 14:22:24 (UTC) on 2016 December 25, near Chiloé, activating the southern part of the Valdivia's earthquake intraplate zone in southern Chile. The epicentre location broadly coincides with the co- and interseismic asperity zone of south of Chiloé Island (Fig. 1 and Moreno *et al.* 2018). The earthquake triggered a small tsunami along the regional littoral (Garetт *et al.* 2018). The Chiloé earthquake process was recorded by multiple instruments, offering an exceptional opportunity to study the first event after the largest earthquake ever recorded. In this
work we developed a fully joint inversion method that incorpo-
rates all available data. This includes data from teleseismic and
strong-motion seismological data, continuous GPS records from
the National Seismological Center of Chile (CSN), tide gauges in-
formation from Hydrographic and Oceanographic Service of the
Chilean Navy (SHOA), and by Sentinel-1 and ALOS-2 radar in-
terferometry. The goal of this research is to better understand the
kinematic source of this event and its relation with the decadal scale
seismic cycle.

## 2 SEISMOTECTONIC SETTING

The austral zone of Chile is located in a complex tectonic con-
text characterized by a subduction margin where two oceanic plates
subduct beneath the South American plate, generating two con-
trasting geodynamic scenarios (Fig. 1). The oceanic plate divergent
contact is controlled by the hot and hydrated Chile Ridge, which,
together with the continental plate, forms the Chile Triple Junction
feature (e.g. Lagabrielle et al. 2000).

The Nazca oceanic plate converges with a gentle obliquity NE-strike (~78°E) with the South American continental plate, at a convergence rate of 6.6 cm yr⁻¹ (Angermann et al. 1999). South of the Taitao Peninsula, the Antarctic oceanic plate also converges obliquely at a SE-strike, with a velocity of 2.1 cm yr⁻¹ (Chase 1978). Oceanic features, such as fracture zones, oceanic plateaux, and ridges may produce local complexities in the subduction plate contact (Lagabrielle et al. 2000; Contreras-Reyes & Carrizo 2011).

The most relevant high oceanic feature in the study zone corre-
sponds to the Guafo Fracture Zone, a transform fault part of the
spreading zone of Chile Ridge that subducts beneath the South
American plate at the south of the Chiloé Island (~43.8°S). This
lithospheric structure strikes N76° E and shows a continuous bathy-
metric scarp that reaches ~1400 km of longitude. The subducted
Nazca plate age decreases from ~14 Ma north of the Guafo Frac-
ture Zone to ~6 Ma to the south, reaching 2 Ma at the Triple-Joint
(Spinelli et al. 2016). Here the subduction-forearc system accom-
modates the oblique convergence by a continental sliver (Chiloé
Block) controlled by the long Liquiñe-Ofqui Fault Zone (LOFZ), a
trench-parallel intra-arc strike-slip fault system of almost 1000 km
long (e.g., Lavenu et al. 1999; Rosenau et al. 2006; Astudillo-
Sotomayor et al. 2021).

## 3 SEISMICITY AND PREVIOUS STUDIES

In the millennium preceding the 1960 Valdivia earthquake, southern Chile was struck by three earthquakes that probably resembled the 1960 main shock (in AD 898–1128, 1300–1398 and 1575) and by five intervening earthquakes, including the events of 1737 and 1837, which ruptured different portions of the plate boundary (Cisternas et al. 2017).

The 1960 Valdivia shock was preceded by a series of foreshocks, which started on 1960 May 21 with an M<sub>w</sub> 8.1 event (Cifuentes 1989; Ojeda et al. 2020), and followed by a large number of after-shocks, most of them located in the northern part of its rupture area, including a slow-slip event (Kanamori & Riviera 2017). Further, the Valdivia event triggered a pronounced drop in the seismicity of the region (Ruiz & Madariaga 2018). The 1960 events were not followed by large events along its rupture area until 2016. Hence, a detailed study of this event provides crucial insights regarding the

way a coupled region will eventually rupture as a mega-earthquake
or smaller event.

Several scientific teams have studied this 2016 event with differ-
ent data and methods. (Fig. 3). Although all the solutions have some
similarities (extent, magnitude), there exist important differences,
such as the maximum slip value (between 2.9 and 5m) or the num-
ber and location of the slip patches. These differences may be due
to several factors, such as the velocity model, the hypocentre loca-
tion, the type of data used, and the smoothing level of the inverse
problem.

For this reason, we have decided to revisit these parameters.
First, we recalculate the hypocentre and moment tensor in order
to get a reliable parametrization of the fault plane. Then, we build
the inverse problem considering all the available data and perform
checkerboard tests to study the resolution of the slip distribution
and the influence of the smoothing level in the solution.

## 4 RELOCATION

Previous studies of the Chiloé earthquake relied on two different hypocentres, one from the USGS, and a second one from the CSN (Table 1). These agencies reported different locations due to the use of different stations, methods, and velocity models. Due to the large discrepancy between these two locations (Inset in Fig. 2), we relocated the epicentre using data from the Chilean National Seismic Network with SeiscompP (GFZ, Gempa) and Hypo71 plugin (Lee et al. 1972). We used the same amount of phases than the CSN solution, and tested two velocity models: the velocity model used by CSN at the time of the event and the 1-D velocity model from Lange et al. (2007). The 1-D velocity model from Lange et al. (2007; Table 2) is the only local model that exists in the region, and it corresponds to the region a further north of the 2016 event. The CSN velocity model is more general for southern Chile (personal communication, Table S1). The relocated hypocentre (Fig. 2) is at 43.50°W (±5 km), 74.21°S (±5 km), and depth of 25 km (±6 km), with an origin time of 14:22:24. The event is southeast of Quilàn Island, closer to the CSN epicentre than to the USGS one.

## 5 MOMENT TENSOR

The moment tensor is determined with the TDMT method (Dreger 2003). Green's functions are calculated using Computer Programs in Seismology (Herrmann 2013). We used the broad-band waveform data from stations of the CSN Network (Supporting Information Fig. S1) at regional distances (from 286 up to 963 km) to estimate the point-source moment tensor. Data are corrected from the instrument response, integrated into displacement, resampled at one sample per second (sp<sub>s</sub>) and bandpass filtered in the frequency band 0.009–0.02 Hz. The horizontal components are rotated to radial and tangential components. Additionally, the source depth was found iteratively by finding the solution that yields the largest variance reduction.

The best solution is shown in Supporting Information Fig. S2 for a centroidal focal depth of 27 km. The moment tensor has 94 per cent of double couples and 6 per cent of compensated linear vector dipole. The variance reduction is 93.3 per cent, which gives a solution of very good quality. The inverted focal mechanism planes have the geometry (strike = 4°, dip = 23°, rake = 99°) for the preferred solution. In addition, we determined the focal mechanism with the FMNEAR method (Delouis 2014) that relies upon local strong-motion data. The preferred plane in this case (0°, 20°, 90°)
<figure><img src="image_4.png" /><figcaption>**Figure 2.** Stations map from strong-motion data (grey), GPS and cGPS (yellow) and tige gauge (red). Epicentral locations: this study (red star), CSN and NEIC (black stars) with a closeup in the inset. Broad-band stations are shown in Supporting Information Fig. S1.</figcaption></figure>

(Supporting Information Fig. S3) is in good accordance with the solution of the TDMT method with a difference of 3° to 4° in the strike and dip. We decided to use the TDMT solution in the inversion because it predicts a larger variance reduction than the FMNEAR mechanism. However, the stations of this region are all located northeast of the event with few stations to the south, which cannot fully constrain the strike. Compared to the NEIC and GCMT solutions, the fault geometry is quite similar, except for the dip value, which is slightly lower (4° and 7° less). The NEIC and GCMT centroids, in contrast to the TDMT centroid, are located to the south of the reported hypocentres. This substantially differs from the centroids of the previous studies (Supporting Information Fig. S4).

# 6 FULLY JOINT INVERSION

## 6.1 Data processing

### 6.1.1 Seismic and GPS data

The Chiloé earthquake was well recorded on global seismological networks as well as by the Chilean National Network of multiparameter stations (broad-band, accelerometer, GPS) and an additional national network of strong-motion data (Fig. 2).

The teleseismic data are the broad-band waveforms of the Federation of Digital Seismograph Networks. We consider teleseismic
<figure><img src="image_5.png" /><figcaption>**Figure 3.** Finite fault solutions from others studies. The used data sets and maximum slip is indicated (SM: strong motion, TEL: Teleseismic body waves. TSU: Mareograph). The yellow star is the epicentre considered in each study. The red star is the centroid of each slip distribution.</figcaption></figure>

<table><caption>**Table 1.** Different locations of the Chiloé earthquake.</caption><thead><tr><td>Agency</td><td>Latitude (°)</td><td>Longitude (°)</td><td>Depth (km)</td></tr></thead><tbody><tr><td>USGS</td><td>-43.416</td><td>-73.951</td><td>35.1</td></tr><tr><td>CSN</td><td>-43.517</td><td>-74.391</td><td>30</td></tr><tr><td>This study</td><td>-43.50</td><td>-74.21</td><td>27</td></tr></tbody></table>

P waves recorded at distances within the range 30°–90° and teleseismic SH waves at distances within the range 30°–75°. After instrument deconvolution, records are integrated to obtain displacement, and then bandpass filtered between 0.01 and 1 Hz for P waves and between 0.01 and 0.4 Hz for SH waves. We use time windows of 80 s around the phase arrival, starting 10 s before P or SH arrival. Teleseismic Green's functions are calculated using generalized ray theory (Helmberger 1983). The strong-motion data are deconvolved from the instrumental response, integrated into velocity, and bandpass filtered in a frequency band 0.02–0.125 Hz. In this case, the Green's function are computed with the spectral discrete wavenumber method of Bouchon (1981), implemented in the Axitra code (Cotton & Coutant 1997). The processed continuous GPS (1 Hz three-component displacement data) are obtained from the global database of high-rate GNSS displacement published by Ruhl *et al.* (2018). The database contains waveforms from the CSN permanent network of high-rate GNSS. We select four stations using signal-to-noise criteria (Fig. 2). Green's function for cGPS are calculated with Axitra code (Cotton & Coutant 1997), and for static GPS with Okada's formulation (Okada 1985).

### 6.1.2 InSAR

We calculate three interferograms from an ascending and a descending track acquired by the C-band Sentinel-1 (S1) satellite and from a descending track of the L-band ALOS-2 satellite (Table 3). The interferograms cover the repeat period of each satellite considering the acquisition observation scenario of the missions, which were 24 and 42 d for S1 and ALOS-2, respectively. We processed the S1 and ALOS-2 interferograms with the ISCE software with standard algorithms for Terrain Observations by Progressive Scans (TOPS) and ScanSAR interferometry (Fattahi *et al.* 2016; Liang & Fielding 2017) respectively. We removed the topographic phase with the 1 arcsec SRTM DEM, and afterwards we averaged pixels 7 times in azimuth and 28 looks in range for the S1 interferograms and 5 looks in range and 2 looks in azimuth for the ALOS-2 interferogram. This mitigates the effect of the temporal decorrelation, increasing the signal-to-noise ratio, resulting in interferograms that have pixels of ~100 m. We removed the ionospheric phase in the ALOS-2 interferogram with the split spectrum algorithm, filtered them with power spectrum coefficients between 0.5 and 0.8, and unwrapped them with the SNAPHU minimum cost flow algorithm (Chen & Zebker 2001). We checked for unwrapping errors in the deforming region in the islands south of the Chiloé island and manually corrected them. Finally, we geocoded the S1 and ALOS-2 interferograms with the 3 arcsec SRTM DEM. Due to the large deformation signal to atmospheric noise ratio, we did not correct the interferograms for atmospheric phase delays. Prior to slip inversion, we removed long-wavelength ramps and we referenced the interferograms to a non-deforming area in the northern part of the Chiloé island. We downsampled the interferograms with a resolution-based algorithm
<table><caption>**Table 2.** Crustal models used for locate and generating synthetics (Lange *et al.* 2007, model)</caption><thead><tr><th>Depth (km)</th><th>V<sub>p</sub> (km s<sup>-1</sup>)</th><th>V<sub>s</sub> (km s<sup>-1</sup>)</th><th>Rho (gr cm<sup>-3</sup>)</th><th>Q<sub>p</sub></th><th>Q<sub>s</sub></th></tr></thead><tbody><tr><td>0</td><td>5.08</td><td>3.14</td><td>2.4</td><td>440</td><td>250</td></tr><tr><td>1</td><td>5.64</td><td>3.17</td><td>2.6</td><td>515</td><td>250</td></tr><tr><td>5</td><td>5.72</td><td>3.33</td><td>2.6</td><td>520</td><td>260</td></tr><tr><td>10</td><td>6.37</td><td>3.8</td><td>2.8</td><td>580</td><td>310</td></tr><tr><td>15</td><td>6.95</td><td>3.87</td><td>3</td><td>680</td><td>320</td></tr><tr><td>23</td><td>7.32</td><td>4.01</td><td>3.1</td><td>760</td><td>340</td></tr><tr><td>30</td><td>7.8</td><td>4.41</td><td>3.2</td><td>900</td><td>430</td></tr><tr><td>39</td><td>7.93</td><td>4.67</td><td>3.3</td><td>940</td><td>500</td></tr><tr><td>45</td><td>7.97</td><td>4.67</td><td>3.3</td><td>960</td><td>500</td></tr><tr><td>60</td><td>8.82</td><td>4.79</td><td>3.5</td><td>1270</td><td>530</td></tr></tbody></table>

<table><caption>**Table 3.** Details of the processed SAR images for interferometry. The columns show the satellite name, radar mode (IW, Interferometric Wide Swath; WD1, ScanSAR), satellite track number, frame numbers, whether the satellite is in an ascending (A) or descending (D) orbit, average incidence angle in the deformation area ($\theta$), radar wavelength ($\lambda$) and date of reference and secondary images (year/month/day).</caption><thead><tr><th>Satellite</th><th>Mode</th><th>Track</th><th>Frames</th><th>Direction</th><th>θ</th><th>λ (cm)</th><th>Reference</th><th>Secondary</th></tr></thead><tbody><tr><td>Sentinel-1</td><td>IW</td><td>83</td><td>730-735-740</td><td>D</td><td>40</td><td>5.55</td><td>2016/12/17</td><td>2017/01/10</td></tr><tr><td>Sentinel-1</td><td>IW</td><td>62</td><td>1034-1039</td><td>A</td><td>35</td><td>5.55</td><td>2016/12/21</td><td>2017/01/14</td></tr><tr><td>ALOS-2</td><td>WD1</td><td>129</td><td>4450-4500</td><td>D</td><td>44</td><td>24.1</td><td>2016/11/26</td><td>2017/01/07</td></tr></tbody></table>

(Lohman & Simons 2005) to increase the density of downsampled points near the deformation signal. This results in ~500–1000 downsampled data points per interferogram which were used in the fault slip inversion.

Green's functions for surface deformation are computed using Okada's formulation (Okada 1985) of the elastic field due to a rectangular dislocation buried in an elastic half-space.

### 6.1.3 Tide gauge

The Chiloé earthquake generated a small tsunami that was recorded by nearby tide stations. We use these data in the slip inversion with appropriate Green's functions calculated with the numerical model NEOWAVE (Yamazaki *et al.* 2009, 2011). NEOWAVE solves the nonlinear shallow water equations in spherical coordinates with non-hydrostatic component in a staggered grid with a shock-capturing finite-difference scheme. In order to save computational resources, NEOWAVE also allows for different nesting grids to reach a sufficiently small grid size to account for inundation processes. The longwave component associated with the tides is removed by subtracting a polynomial fitting, and high-frequency noise is suppressed with a three-point rolling average (Adriano *et al.* 2018). The resulting waveforms define the target tsunami data to be inverted. Next, for each element of the fault discretization, NEOWAVE computes the static seafloor displacement and uses it as the tsunami initial condition. This is based on Okada's equations (Okada 1985) and the horizontal contributions (Tanioka & Satake 1996) of a unit slip source. Since tide gauges are sampled each 60 s, and due to the relatively short duration of the earthquake rupture, source kinematics can be safely neglected in the tsunami modelling. The outputs from the Green functions are only used in the first peak-to-peak portion of the observed signal to keep the superposition principle, avoiding the nonlinear effects of tsunami waves with the shore (Williamson *et al.* 2018).

The particular case of the Castro station requires detailed bathymetry data to correctly simulate the tsunami propagation inside the Castro fjord. This is achieved with a three-level nesting. Bathymetry data with a resolution of 15 arcsec is from the General

Bathymetry Chart of the Ocean (GEBCO), while high-resolution DEMs are built with the global 1 arcsec topography data from the USGS and the nautical charts provided by the Servicio Hidrográfico y Oceanográfico de la Armada (SHOA, Chilean Navy).

The Level 1 (coarsest) data are gridded with a resolution of 30 arcsec, resampling the original 15 arcsec data. Level 2 is nested in level 1 with 10 arcsec, and finally, Level 3 is nested in level 2 with 2 arcsec. Both levels 2 and 3 are interpolated from the 15 arcsec bathymetry, 1 arcsec topography and the nautical charts (Supporting Information Fig. S5). Finally, each tsunami simulation is computed with 5 hr of propagation time, summing up $204 \times 5 = 1020$ Green's function.

## 6.2 Methodology of the joint inversion

We apply a linear waveform inversion scheme with a finite-fault model. A grid of point sources is distributed uniformly over the fault plane, with dimensions substantially larger than the expected extent of the rupture. The point sources are triggered according to the passage of a circular rupture front propagating at a prescribed speed. We adopt the time windows method (Olson & Apsel 1982; Hartzell *et al.* 1983), in which the slip-rate function on each subfault is expanded into a series of overlapping isosceles triangle functions with equal rise time. In the case of joint inversion using independent slip-time history data (GPS, Tsunami and InSAR) and dependent slip-time history data (body waves, cGPS and strong motion), this temporal rupture model is preserved for the static case for consistency. For an inversion that relies upon independent slip-time history data only, we use one isosceles triangle to recover the coseismic slip of each point source.

Green's functions are computed for two slip directions, so the final slip of each subfault will be a linear combination of these directions. The matrix of synthetics **G** is constructed placing end-to-end synthetic Green's functions for all stations for each subfault for a triangle source-time function of a given rise time. **d** is the concatenated observed data at all stations. The inversion problem is expressed as the least-squares problem $Gm = d$, in which **m** is a vector containing the unknown moment of each subfault. In
addition, a priori data covariance matrix C is used to scale the system equations. Appending a smoothing operator L to the linear system, the inverse problem is expressed as:

$$ \begin{bmatrix} C^{-1}G \\ \lambda L \end{bmatrix} \mathbf{m} = \begin{bmatrix} C^{-1}d \\ 0 \end{bmatrix}, \qquad (1) $$

in which $\lambda$ is a regularization weighting parameter. The smoothing parameter value is selected as a trade-off between minimizing the model misfit and minimizing the solution roughness. We impose the no backslip condition so $m_i \ge 0$, for all i (for each subfault). We find the best model m by non-negative least-squares methods (Lawson & Hanson 1974). Given the different resolution levels of the data sets, different values of the $\lambda$ parameter are necessary for each individual inversion (Supporting Information Fig. S6). The inverse problem for two different data sets is represented by the following two linear systems

$$ \begin{bmatrix} G_1 \\ \lambda_1 L \end{bmatrix} \mathbf{m} = \begin{bmatrix} d_1 \\ 0 \end{bmatrix} \quad \text{and} \quad \begin{bmatrix} G_2 \\ \lambda_2 L \end{bmatrix} \mathbf{m} = \begin{bmatrix} d_2 \\ 0 \end{bmatrix}. \qquad (2) $$

Following Kaverina et al. (2002), we multiply the first system equation by $\lambda_2$ and the second one by $\lambda_1$. The two systems may be combined into a single system of equations

$$ \begin{bmatrix} \lambda_2 G_1 \\ \lambda_1 G_2 \\ \lambda_1 \lambda_2 L \end{bmatrix} \mathbf{m} = \begin{bmatrix} \lambda_2 d_1 \\ \lambda_1 d_2 \\ 0 \end{bmatrix}. \qquad (3) $$

We repeat this procedure for all data sets. The smoothing parameter of the joint inversion is obtained as $\lambda_{joint} = \prod_i \lambda_i$, and the weight of ith data set is $w_i = \prod_{j \neq i} \lambda_j$.

# 7 RESULTS

## 7.1 Checkerboard tests

In order to evaluate the resolution of the inversion provided by different data sets, a checkerboard test is performed. We consider a fault plane with geometry (0,20,99) to construct the synthetic slip distribution, with a rupture velocity of 2.8 km s⁻¹. The synthetic data are generated for the same stations that provided the real data, and are then inverted using the same scaling factor and smooth parameters determined to the data sets in the main-shock inversion. Finally, we perform the inversion using a rupture velocity of 3 km s⁻¹ over a fault plane with geometry (4,23,99) to simulate some error.

To test the overall resolution of each data set over the whole fault plane, two checkerboard test with patch size of 60 km are performed (Fig. 4). In both tests, teleseismic, cGPS, GPS and strong motion, InSAR and tide gauge data have limited resolution over the whole plane. The joint inversion is the one that best reproduces the input fault slip from the fault updip to downdip and along the fault strike.

Fig. 5 shows the sensibility test of the slip distribution with respect to the smoothing parameter. The slip patches form is not mainly sensitive to this parameter. The main differences are found in the maximum slip obtained in the inversions. A smoothing parameter lower than $\lambda_{joint}$ generates an overestimation of maximum slip in all cases tested. Increasing the smoothing biases the maximum slip in the overestimated models, but could slightly affect the overall shape of the slip distribution.

To estimate the minimum size of patches recovered by the joint inversion, checkerboard tests with patch size of 40 and 30 km are performed (Fig. 6). The joint inversion has high sensitivity to the

updip slip distribution and to the slip distribution near the hypocentre. In this zone, it is possible to recover the slip patches of 30 km. The slip in the downdip region is either underestimated or slightly displaced.

## 7.2 Coseismic slip distribution

Teleseismic, strong motion, GPS, cGPS, InSAR and tide gauge data are jointly inverted for the 2016 December 25 event to retrieve the coseismic slip distribution. The fault has a strike of 4°, a dip of 23° and a variable rake. The size of the fault is 170 × 120 km, divided into 10 × 10 km subfaults. The result of the inversion is shown in Fig. 7. The final slip is composed of a relatively simple slip patch located south of Chiloé Island with a maximum slip of 3.2 m close to the hypocentre. The maximum rupture velocity is 3.2 km s⁻¹. In comparison, the slip of each individual data sets inversions are shown in Supporting Information (Fig. S7). Each individual slip model shows a smoother pattern compared to that of the joint inversion. The main characteristics of these models are reported in Supporting Information Table S2. The rupture propagates mainly updip over the seismogenic zone with an extension in both northeast and southwest directions. Except for the small strip of slip going to the south, the rupture seems not to propagate further south than the Guafo Island. This slip distribution produces a good overall fit for each seismic, geodetic, and tsunami data (Figs 8–11).

We perform several inversions to search for a range of maximum slip. The maximum amount of slip which explains the data is obtained by solving the inverse problem without regularization ($\lambda = 0$). The unregularized slip distribution is shown in Fig. 12(a), and shows the maximum slip patch with 4 m closer to the hypocentre. The slip distribution is similar to the regularized solution Fig. 12(b), except for the slip patch of the maximum slip. In this manner, the $\lambda$ parameter in the range [0, $\lambda_{joint}$] controls the amount of slip without degrading the rupture pattern (Fig. 12d). To explore lower values of maximum slip without degrading the data fit, we perform different joint inversions for $\lambda = \lambda_{joint}$ adding the restriction $0 \le m_i \le \max_{\text{slip}}$ for different values of maximum slip. We find solutions with maximum slip between 3.2 and 2.5 m that fits the data equally well (Supporting Information Fig. S9). Solution with maximum slip of 2.5 m is shown in Fig. 12(c). This results in a range of maximum slip of 3.2 ± 0.8 m.

# 8 DISCUSSION

## 8.1 Joint inversion

The inversion resolution is linked to the formulation of the problem, and on the data set. The main difficulty in a joint inversion is how to determine the relative weights for each data set in such a way that only the components of each data that have the resolution to resolve a specific aspect of the rupture process are incorporated in the slip model. In general, only certain characteristics such as the location and timing of large slip, average direction and rupture velocity are well resolved because fault slip models have only limited reliability due to the implicit underdetermined problem. To obtain a reliable image of the details of the source process, many conditions are required: high spatial and temporal resolution data with a good spatial coverage, records of good quality, a good parametrization of the source model, and good knowledge of the velocity model. Unfortunately, the study area in southern Chile is characterized by large fjords, resulting in difficult field logistics, and the distribution
Checkerboard Test 1

<figure><img src="image_6.png" /><figcaption>**Figure 4.** Checkerboard tests for the inversions. Slip distribution obtained using the different data sets: (a) input model, (b) joint inversion, (c) teleseismic, (d) strong motion, (e) GPS, (f) cGPS, (g) InSAR and (h) Tide gauge. Slip scale ranges from blue (0 m) to red (3.2 m). The star is the epicentre.</figcaption></figure>

of seismic and geodetic station is very sparse. Hence, neither the geodetic nor the strong-motion stations uniformly cover the near-field region, limiting the resolution for an inversion based on each individual data set. To overcome this limitation, we have built a joint inversion to take advantage of each data set. Each one brings a different piece of information about the spatial, temporal evolution, and moment of the source (Yue et al. 2020). Seismic waveform data are essential to resolve the temporal evolution of fault slip.

Checkerboard Test 2

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

Strong-motion records are especially important to obtain a fine image of rupture propagation for large earthquakes, but need a very good station distribution. For seismic data, the main problem are the uncertainties that propagate from the inaccuracy of Green's functions due to an uncertain velocity model. Even though nearfield seismic data are essential to resolve the temporal change in detail, a trade-off between slip location and rupture time exists in inversion when using only waveform data (Koketsu 2016). Geodetic data are
<figure><img src="image_8.png" /><figcaption>**Figure 5.** Sensibility of Checkerboard tests with respect to smoothing parameter.</figcaption></figure>

important for slip inversion to resolve the final static slip distribution (e.g. Ward & Barrientos 1986; Yabuki & Matsu'ura 1992). Geodetic data are combined with seismic data to better constrain the amount of total slip, which is often poorly constrained from seismic data alone. Finally, tide gauge records of tsunamis excited by

earthquakes provide information about the low-frequency behaviour of the spatial distribution of fault slip. In the case of the inversion of the Chiloé event, the checkerboard tests show that the fully joint inversion improves significantly the resolution of the slip distribution form. Even the solution obtained without regularization does not
<figure><img src="image_9.png" /></figure>

<figure><img src="image_10.png" /><figcaption>**Figure 6.** Checkerboard inversions for patch size of 40 and 30 km.</figcaption></figure>

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

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

<figure><img src="image_13.png" /><figcaption>**Figure 7.** Joint inversion final slip distribution along strike and dip, showing the discretization of the fault (10 km × 10 km subfaults). The maximum slip is 3.2 m north west of the epicentre (blue star). The Source Time Function (STF) is shown in the inset. The total rupture time is 40 s. Black arrows show the variable rake directions.</figcaption></figure>

<figure><img src="image_14.png" /><figcaption>**Figure 8.** Horizontal and vertical coseismic static displacement field at GPS sites. Observed (black/white arrows) and synthetic (red arrows) for the stations shown in Fig. 2. Contour lines of the final slip distribution are drawn from blue to red.</figcaption></figure>

present spurious slip updip, as expected on finite-fault inversions of megathrust earthquakes offshore. The checkerboard tests indicate that the resolving power of our data is high in the middle and upper part of the fault plane, where the slip distribution of the Chiloé event is located. In this area, the minimum size of the slip patches that can be recovered in the inversion is nearly 30 km. However, the fully joint inversion is not enabled to resolve the maximum slip. We found slip distributions with maximum slip between 4 and 2.5 m which explain the data with a similar fitting level.
<figure><img src="image_15.png" /><figcaption>**Figure 9.** Far field body wave fit for the joint inversion. Observed (black lines) and synthetic (red lines) for *P* and *SH* displacements waveforms. The azimuth and distance in degrees are indicated at the beginning of each record with the azimuth on top.</figcaption></figure>

<figure><img src="image_16.png" /><figcaption>**Figure 10.** Near field waveform fit for the stations shown in Fig. 2. (Observed are black and model predictions are red). (a) Strong motion. Number below of each trace is the peak velocity. (b) Continuous GPS. (c) Tide gauge. The shaded area in the time-series depict the trace of the inverted signal.</figcaption></figure>

## 8.2 Slip deficit in the Valdivia rupture segment

Recent observations indicate that the seismic coupling patterns may change in time (e.g. Bedford *et al.* 2020; Mallick *et al.* 2021). Hence, the model of slip budget based on the balance of energy accumulation as slip deficit and energy release as coseismic slip is in general too simplistic. This is especially problematic in areas that have had a lack of continuous observations like the Chiloé segment, where a sparse network of continuous GPS data has been available only since 2010–2015. As such, these transient changes might have been completely unnoticed. For example, Luo *et al.* (2020) have shown

that the plate coupling pattern 250 km south of the study area changed a few years prior to 2010 due to an hypothetical slow-slip event based on the analysis of a single GNSS station. We acknowledge that the aforementioned model of slip budget is a simplistic approximation of the earthquake cycle, but the lack of long-term time-series of strain accumulation in the study area precludes a more thorough analysis. In the following we only use the maximum slip in the study area from Moreno *et al.* (2009) who inverted the costal uplift/subsidence data of Plafker & Savage (1970), and neglect sources of uncertainty due to aseismic deformation because they are not known.
<figure><img src="image_17.png" /><figcaption>**Figure 11.** InSAR maps for the joint inversion: (a) Downsampled observed deformation, (b) model prediction and (c) residuals comparing data and the model for the inverted slip. Models are computed for three interferograms from an ascending and a descending track acquired by the C-band Sentinel-1 (S1) satellite and from a descending track of the L-band ALOS-2 satellite.</figcaption></figure>

Previous studies have addressed the earthquake potential in the Valdivia segment comparing the coseismic slip of the 1960 and 2016 earthquakes with the estimated slip deficit accumulated since 1960. The convergence of the Nazca plate has a geodetic rate of 6.6 cm yr⁻¹ (Angermann et al. 1999), resulting in a slip deficit of ~3.7 m in the 56 years between the two earthquakes. Moreno et al. (2011) showed that the megathrust was almost completely locked in the ruptured area, albeit the locking models have limited along strike spatial resolution due to the lack of a dense GNSS monitoring network in the 2016 rupture area. Previous estimates of coseismic peak slip of the 2016 earthquake are 4.5 m (Ruiz et al. 2017), 4.7 m (Xu 2017), 5 m (Melgar et al. 2017), 4.2 m (Lange et al. 2018) and 2.9 m (Moreno et al. 2018). All these slip estimates, except for the latter, are larger than the accumulated slip deficit since the 1960 rupture, suggesting that the Chiloé earthquake released some slip that was accumulated before the 1960 earthquake (Melgar et al. 2017; Xu 2017; Lange et al. 2018). This led Xu (2017) to suggest that the Chiloé segment also behaves like that of the 2016 Pedernales earthquake in Ecuador where earthquake supercycles release more slip than the slip accumulated since the previous rupture (Nocquet et al. 2017). The previous earthquake of equivalent magnitude in the Chiloé segment occurred in 1575 (Cisternas et al. 2005). Following the same reasoning, the 1575–1960 interval accumulated ~25 m of slip deficit in an asperity that ruptured 20–25 m in the Valdivia earthquake (Moreno et al., 2009) under the assumption of no transient changes in plate coupling. Hence, the 1960 event should have released most of the energy accumulated between that earthquake and the previous event in 1575. Our slip distribution provides insights in this aspect. Our slip model predicts a similar

slip patch in magnitude and location to that of Moreno *et al.* (2018) who used the same GNSS and InSAR data that we used, but without seismological data. The inversion smoothing introduces uncertainties of ~0.8 m in our slip model. Considering these uncertainties the maximum magnitude of our slip model overlaps with that of Moreno *et al.* (2018). Our joint slip model predicts a maximum slip of ~3.2 ± 0.8 m, differ by 9 per cent with respect to that of Moreno et al., (2018), and are at least 25 per cent lower than the other estimates (Melgar *et al.* 2017; Ruiz *et al.* 2017; Xu 2017; Lange *et al.* 2018). Hence, maximum slip in the 2016 asperity is in general not larger than the slip deficit accumulated since 1960. Therefore, under the assumption of no occurrence of a slow-slip earthquake, it is unlikely that the Chiloé earthquake released strain accumulated before the 1960 earthquake in the 2016 rupture area. In contrast, other studies (Melgar *et al.* 2017; Ruiz *et al.* 2017; Xu 2017; Lange *et al.* 2018) achieved the conclusion that the Chiloé earthquake released energy accumulated before 1960 by not considering a joint inversion where all the different data sets were properly tested for their resolving power, and taking into account regularization uncertainties.

A comparison between our slip model with a high resolving power with the locking distribution of Moreno *et al.* (2011) is not direct because the latter has limited resolution immediately south of Guafo island and in the offshore zones west of the Chiloé Island. Hence, the rupture model cannot explain why the rupture stopped near Guafo Island and why it did not propagate to the north of the Chiloé Island. The latter is an area with a highly coupled patch which should have acted as an asperity that promoted the rupture propagation instead of a barrier (Moreno *et al.* 2018).
<figure><img src="image_18.png" /><figcaption>**Figure 12.** (a) Slip distribution solution without regularization constrain. (b) Preferred finite-fault solution. (c) Slip distribution solution adding the restriction $0 \le m_i \le 2.5$ m. (d) Comparison of the maximum slip along strike, for the slip solutions above.</figcaption></figure>

## 9 CONCLUSION

The joint inversion of teleseismic, regional seismological strong motion, static and continuous GNSS, InSAR and tide gauge data results in coseismic slip models that have a high resolving power. Hence, the slip model have a high confidence that the maximum coseismic slip of the Chiloé earthquake was located next to the southwest corner of the island, in the downdip extent of the megathrust, with a maximum slip of $3.2 \pm 0.8$ m. The previous suggestions that the Chiloé earthquake released seismic moment accumulated before the Valdivia earthquake are not supported in the study area by the integration of the multiparametric data set. This implies that even in areas with relatively sparse instrumentation like southern Chile it is possible to overcome the model resolution limitations of each individual data set by incorporating a diverse multiparametric data set, resulting in robust estimates of fault slip. We speculate that previous studies might be biased by considering limited data sets, and that fault slip models of earthquakes that occur in areas with poor data coverage like southern Chile will have to rely on the complete ensemble of seismological and geodetic earthquake data. Although our inversion benefited from the fact that the rupture area was located above areas with geodetic measurements, we speculate that a multiparametric joint inversion can also improve fault slip estimates for ruptures in the trench, in which geodetic data have low resolving power. Finally, another aspect that is important to highlight is why the first rupture after the Valdivia earthquake occurred in a deep section of the megathrust instead of a section closer to the trench. Our slip model provides no insights in this aspect, but frictional models should be explored in the future, in particular the

role of the pore fluid pressure, and how it modulates the interseismic strain build-up (e.g. Moreno *et al.* 2018). Future studies should also address what arrested the earthquake rupture propagation towards its southern and northern parts, areas that were highly coupled but that did not slip during the earthquake.

## ACKNOWLEDGMENTS

We thank the Servicio Hidrográfico y Oceanográfico de la Armada de Chile (SHOA) for providing the nautical charts 7372 and 7370, IOC for tide gauges records, and the CSN for geodetic and seismological data. Also, we are grateful to the international networks that provide the rest of the data accessible through the Incorporated Research Institutions for Seismology (IRIS). Moment tensors were computed using the mtpackagev3.1 package developed by Douglas Dreger and Sean Ford of the Berkeley Seismological Laboratory, and Green's functions were computed using the Computer Programs for Seismology software developed by Robert Herrmann. We thank B. Delouis for the FMNEAR code. We thank D. Melgar, M. Moreno, D. Lange, S. Ruiz and W. Xu for providing their finite fault solution, and R. Grandin for discussions on joint inversions. This work was supported by the PRS project (https://prs.dgf.uchile.cl/). F.D. was funded by FONDECYT de Iniciación 11220513 research grant.

## DATA AVAILABILITY

Records from Federation of Digital Seismograph Networks can be obtained from IRIS Data Management Center at www.iris.edu.
Strong-motion records were obtained from http://evtdb.csn.uchile.cl, which is maintained by Centro Sismológico Nacional (CSN). The high-rate GNSS displacement waveforms belong to the data set published in Seismological Research Letters by Ruhl *et al.* (2018), available at http://at.dx.doi.org/10.5281/zenodo.1434374. Sentinel-1 data are property of the European Space Agency and are open-source data available for download at Vertex by the Alaska Satellite Facility and Plateforme d'Exploitation des Produits SentinEL (PEPS) by CNES. ALOS-2 data were provided by JAXA through the scientific quota of the CEOS Latin America Volcano Pilot Project coordinated by Michael Poland (U.S. Geological Survey). The ALOS-2 SLC data are available in the UNAVCO Seamless SAR Archive SSARA (https://web-services.unavco.org/brokered/ssara/gui). Tide gauges records are downloaded from the IOC website (http://www.ioc-sealevelmonitoring.org/). The recorded displacements at GUAf and MELK stations are from Moreno *et al.* (2018). 15 arcsec bathymetry data were downloaded from GEBCO website (https://www.gebco.net/data_and_products/gridddęd_bathymetry_data/).

## REFERENCES

Adriano, B., Hayashi, S. & Koshimura, S., 2018. Analysis of spatio-temporal tsunami source models for reproducing tsunami inundation features, *Geosciences*, **8**(1), 3, doi:10.3390/geosciences8010003.

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.

Astudillo-Sotomayor, L., Jara-Muñoz, J., Melnick, D., Cortés-Aranda, J., Tassara, A. & Strecker, M. R., 2021. Fast Holocene slip and localized strain along the Liquiñe-Ofqui strike-slip fault system, Chile, *Sci. Rep.*, **11**, 5970, doi:10.1038/s41598-021-85036-5.

Barrientos, S. & Ward, S., 1990. The 1960 Chile earthquake: inversion for slip distribution from surface deformation, *Geophys. J. Int.*, **103**, 589-598.

Bedford, J.R., Moreno, M., Deng, Z., Oncken, O., Schurr, B., John, T., Báez, J.C. & Bevis, M., 2020. Months-long thousand-kilometre-scale wobbling before great subduction earthquakes, *Nature*, **580**, 628-635.

Bouchon, M., 1981. A simple method to calculate Green's functions for elastic layered media, *Bull. seism. Soc. Am.*, **71**, 959-971.

Chase, C., 1978. Plate kinematics: the Americas, East Africa, and the rest of the World, *Earth planet. Sci. Lett.*, **37**(3), 355-368.

Chen, C. W. & Zebker, H. A., 2001. Two-dimensional phase unwrapping with use of statistical models for cost functions in nonlinear optimization, *J. Opt. Soc. Am. A*, **18**(2), 338-351.

Cifuentes, I., 1989. The 1960 Chilean earthquakes, *J. geophys. Res.*, **94**(B1), 665-690.

Cisternas, M. *et al.*, 2005. Predecessors of the giant 1960 Chile earthquake, *Nature*, **437**(7057), 404-407.

Cisternas, M., Garrett, E., Wesson, R., Dura, T. & Ely, L., 2017. Unusual geologic evidence of coeval seismic shaking and tsunamis shows variability in earthquake size and recurrence in the area of the giant 1960 Chile earthquake, *Mar. Geol.*, **385**, 101-113.

Contreras-Reyes, E. & Carrizo, D., 2011. Control of high oceanic features and subduction channel on earthquake ruptures along the Chile-Peru subduction zone, *Phys. Earth planet. Inter.*, **186**, 49-58.

Cotton, F. & Coutant, O., 1997. Dynamic stress variations due to shear faults in a plane-layered medium, *Geophys. J. Int.*, **128**, 676-688.

Delouis, B., 2014. FMNEAR: determination of focal mechanism and first estimate of rupture directivity using near-source record and a linear distribution of point sources, *Bull. seism. Soc. Am.*, **104**(3), 1479-1500.

Dreger, D. S., 2003. TDMT_INV: time domain seismic moment tensor inversion, in *International Handbook of Earthquake and Engineering Seismology*, Vol. 81B, pp. 1627, Academic Press.

Duputel, Z. *et al.*, 2015. The Iquique earthquake sequence of April 2014: Bayesian modeling accounting for prediction uncertainty, *Geophys. Res. Lett.*, **42**, doi:10.1002/2015GL065402.

Fattahi, H., Agram, P. & Simons, M., 2016. A network-based enhanced spectral diversity approach for TOPS time-series analysis, *IEEE Trans. Geosci. Remote Sens.*, **55**, 777-786.

Garrett, E., Brader, M., Melnick, D., Bedford, J. & Aedo, D., 2018. First field evidence of coseismic land-level change associated with the 25 December 2016 Mw 7.6 Chiloe, Chile, earthquake, *Bull. seism. Soc. Am.*, **109**(1), 87-98.

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**(6), 1553-1583.

Helmberger, D. V., 1983. Theory and application of synthetic seismograms in earthquakes: observation, theory and interpretation, in *Proc. Int. Sch. Phys 'Enrico Fermi', Course LXXX*, eds Kanamori, H. & Boshi, E., pp. 174-221, North Holland Publishing Company.

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

Kanamori, H., 1977. The Energy Release in Great Earthquakes, *J. geophys. Res.*, **82**(20), 2981-2987.

Kanamori, H. & Riviera, L., 2017. An Mw=7.7 slow earthquake in 1960 near the Aysén Fjord region, Chile, *Geophys. J. Int.*, **211**, 93–106.

Kaverina, A., Dreger, D. & Price, E., 2002. The combined inversion of seismic and geodetic data for the source process of the 16 October 1999 Mw 7.1 Hector Mine, California, earthquake, *Bull. seism. Soc. Am.*, **92**, 1266-1280.

Koketsu, K., 2016. An overview of joint inversion in earthquake source imaging, *J. Seismol.*, **20**(4), 1131-1150.

Lagabrielle, Y., Guivel, C., Maury, R. C., Bourgois, J., Fourcade, S. & Martin, H., 2000. Magmatic-tectonic effects of high thermal regime at the site of active ridge subduction: the Chile Triple Junction model, *Tectonophysics*, **326**, 255-268.

Lange, D., Rietbrock, A., Haberland, C., Bataille, K., Dahm, T., Tilmann, F. & Flüh, E. R., 2007. Seismicity and geometry of the south Chilean subduction zone (41.5°S-43.5°S): Implications for controlling parameters, *Geophys. Res. Lett.*, **34**(6), L06311.1-L06311.5.

Lange, D., Ruiz, J., Carrasco, S. & Manríquez, P., 2018. The Chiloe Mw 7.6 earthquake of 2016 December 25 in Southern Chile and its relation to the Mw 9.5 1960 Valdivia earthquake, *Geophys. J. Int.*, **213**, 210-221.

Lavenu, A. & Cembrano, J., 1999. Estados de esfuerzo compresivo plioceno y compresivo-transpresivo pleistoceno, Andes del sur, Chile (38-42°30'S), *Revista geológica de Chile*, **26**(1), 67-87.

Lawson, C. L. & Hanson, B. J., 1974. *Solving Least Squares Problems*, Prentice-Hall.

Lee, W. H. K. & Lahr, J. C., 1972. HYPO71: a computer program for determining hypocenter, magnitude, and first motion pattern of local earthquakes, Open File Report 72-224.

Liang, C. & Fielding, E.J., 2017. Measuring azimuth deformation with L-band ALOS-2 ScanSAR interferometry, *IEEE Trans. Geosci. Remote Sens.*, **55**(5), 2725-2738.

Lohman, R. B. & Simons, M., 2005. Some thoughts on the use of InSAR data to constrain models of surface deformation: noise structure and data downsampling, *Geochem. Geophys. Geosyst.*, **6**(1), Q01007, doi:10.1029/2004GC000841.

Luo, H., Ambrosius, B., Russo, R. M., Mocanu, V., Wang, K., Bevis, M. & Fernandes, R., 2020. A recent increase in megathrust locking in the southernmost rupture area of the giant 1960 Chile earthquake, *Earth planet. Sci. Lett.*, **537**, 116200.

Maldonado, V., Contreras, M. & Melnick, D., 2021. A comprehensive database of active and potentially-active continental faults in Chile at 1:25,000 scale, *Sci. Data*, **8**, 20, doi:10.1038/s41597-021-00802-4.

Mallick, R., Meltzner, A. J., Tsang, L. L., Lindsey, E. O., Feng, L. & Hill, E. M., 2021. Long-lived shallow slow-slip events on the Sunda megathrust, *Nat. Geosci.*, **14**(5), 327-333.

Meade, B. J., Klinger, Y. & Hetland, E. A., 2013. Inference of multiple earthquake-cycle relaxation timescales from irregular geodetic sampling of interseismic deformation, *Bull. seism. Soc. Am.*, **103**(5), 2824-2835.
Melgar, D., Riquelme, S., Xu, X., Baez, J. C., Geng, J. & Moreno, M., 2017. The first since 1960: a large event in the Valdivia segment of the Chilean Subduction Zone, the 2016 M7.6 Melinka earthquake, *Earth planet. Sci. Lett.*, **474**, 68–75.

Melnick, D., et al., 2018. Back to full interseismic plate locking decades after the giant 1960 Chile earthquake, *Nat. Commun.*, **9**, 3527, doi:10.1038/s41467-018-05989-6.

Moreno, M. S., Bolte, J., Klotz, J. & Melnick, D., 2009. Impact of megathrust geometry on inversion of coseismic slip from geodetic data: application to the 1960 Chile earthquake, *Geophys. Res. Lett.*, **36**, L16310, doi:10.1029/2009GL039276.

Moreno, M. et al., 2011. Heterogeneous plate locking in the South-Central Chile subduction zone: building up the next great earthquake, *Earth planet. Sci. Lett.*, **305**, 413–424.

Moreno, M. et al., 2018. Chilean megathrust earthquake recurrence linked to frictional contrast at depth, *Nat. Geosci.*, **11**, 285–290.

Nocquet, J.M. et al., 2017. Supercycle at the Ecuadorian subduction zone revealed after the 2016 Pedernales earthquake, *Nat. Geosci.*, **10**, 145–149.

Ojeda, J., Ruiz, S., del Campo, F. & Carvajal, M., 2020. The 21 May 1960 Mw 8.1 Concepción earthquake: a deep megathrust foreshock that started the 1960 central–south Chilean seismic sequence, *Seismol. Res. Lett.*, **91**(3), 1617–1627.

Okada, Y., 1985. Surface deformation due to shear and tensile faults in a half-space, *Bull. seism. Soc. Am.*, **75**(4), 1135–1154.

Olson, A. H. & Apsel, R. J., 1982. Finite faults and inverse theory with applications to the 1979 Imperial Valley earthquake, *Bull. seism. Soc. Am.*, **72**(6A), 1969–2001.

Plafker, G. & Savage, J., 1970. Mechanism of the Chilean earthquakes of May 21 and 22, 1960, *Bull. geol. Soc. Am.*, **81**(4), doi:10.1130/0016-7606(1970)81[1001:MOTCEO]2.0.CO;2.

Rosenau, M., Melnick, D. & Echter, H., 2006. Kinematic constraints on intra-arc shear and strain partitioning in the southern Andes between 38°S and 42°S latitude, *Tectonics*, **25**(4), doi:10.1029/2005TC001943.

Ruhl, C. J., Melgar, D., Geng, J., Goldberg, D. E., Crowell, B. W., Allen, R. M. & d'Anastasio, E., 2019. A global database of strong-motion displacement GNSS recordings and an example application to PGD scaling, *Seismol. Res. Lett.*, **90**(1), 271–279.

Ruiz, S., Moreno, M., Melnick, D., del Campo, F., Poli, P., Baez, J. C., Leyton, F. & Madariaga, R., 2017. Reawakening of large earthquakes in south central Chile: The 2016 Mw 7.6 Chiloé event, *Geophys. Res. Lett.*, **44**, 6633–6640.

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

Satake, K. & Kanamori, H., 1991. Use of tsunami waveforms for earthquake source study, *Nat. Hazards*, **4**(2–3), 193–208.

Spinelli, G. A., Wada, I., He, J. & Perry, M., 2016. The thermal effect of fluid circulation in the subducting crust on slab melting in the Chile subduction zone, *Earth planet. Sci. Lett.*, **434**, 101–111.

Tanioka, Y. & Satake, K., 1996. Tsunami generation by horizontal displacement of ocean bottom, *Geophys. Res. Lett.*, **23**(8), 861–864.

Xu, W., 2017. Finite-fault slip model of the 2016 Mw 7.5 Chiloé earthquake, southern Chile, estimated from Sentinel-1 data, *Geophys. Res. Lett.*, **44**, 4774–4780.

Yabuki, T. & Mats'u'ra, M., 1992. Geodetic data inversion using a Bayesian information criterion for spatial distribution of fault slip, *Geophys. J. Int.*, **109**(2)(2), 363–375.

Yamazaki, Y., Kowalik, Z. & Cheung, K. F., 2009. Depth-integrated non-hydrostatic model for wave breaking and run-up, *Int. J. Numer. Methods Fluids*, **61**(5), 473–497.

Yamazaki, Y., Cheung, K. F. & Kowalik, Z., 2011. Depth-integrated, non-hydrostatic model with grid nesting for tsunami generation, propagation, and run-up, *Int. J. Numer. Methods Fluids*, **67**(12), 2081–2107.

Yokota, Y., Koketsu, K., Fujii, Y., Satake, K., Sakai, S., Shinohara, M. & Kanazawa, T., 2011. Joint inversion of strong motion, teleseismic, geodetic, and tsunami datasets for the rupture process of the 2011 Tohoku earthquake, *Geophys. Res. Lett.*, **38**(7), doi:10.1029/2011GL050098.

Yoshida, S. & Koketsu, K., 1990. Simultaneous inversion of wave-form and geodetic data for the rupture process of the 1984 Naganoken-Seibu, Japan, earthquake, *Geophys. J. Int.*, **103**, 355–362.

Yue, H., Lay, T., Rivera, L., An, C., Vigny, C., Tong, X. & Baez Soto, J. C., 2014. Localized fault slip to the trench in the 2010 Maule, Chile Mw = 8.8 earthquake from joint inversion of high-rate GPS, teleseismic body waves, InSAR, campaign GPS, and tsunami observations, *J. geophys. Res.*, **119**, 7786–7804.

Yue, H., Zhang, Y., Gai, Z., Wang, T. & Zhao, L., 2020. Resolving rupture processes of great earthquakes: reviews and perspective from fast response to joint inversion, *Sci. China Earth Sci.*, **63**(4), 492–511.

Wald, D. J. & Heaton, T. H., 1994. Spatial and temporal distribution of slip for the 1992 Landers, California, earthquake, *Bull. seism. Soc. Am.*, **84**, 668–691.

Ward, S. N. & Barrientos, S. E., 1986. An inversion for slip distribution and fault shape from geodetic observations of the 1983, Borah Peak, Idaho, Earthquake, *J. geophys. Res.*, **91**(B5), 4909, doi:10.1029/jb091ib05p04909.

Williamson, A. L. & Newman, A. V., 2018. Limitations of the resolvability of finite-fault models using static land-based geodesy and open-ocean tsunami waveforms, *J. geophysical Res.*, **123**(10), 9033–9048.

## SUPPORTING INFORMATION

Supplementary data are available at GJI online.

**Figure S1.** Map of the broad-band stations from CSN and Geoscope used in the moment tensor inversion. The focal mechanism found is plotted at the epicentre location. Geoscope station was also used in FMNEAR.

**Figure S2.** Moment tensor inversion. The data are shown as black lines and the synthetics are shown as red lines. For each station, three-component data and synthetics are compared, and the distance, the azimuth, the maximum three-component trace amplitude and fitting parameter such as the variance and the variance reduction (Var. Red.) are provided. Solution information includes the strike, rake and dip for the two possible double-couple (DC) planes, the scalar seismic moment and $M_w$. Information about the moment tensor decomposition in terms of percent DC, compensated linear vector dipole (CLVD) and isotropic (ISO) are also listed.

**Figure S3.** Results obtained with the real-time FMNEAR (for determination of focal mechanism using NEAR-source records) algorithm implemented in Chile. In the case of the 2016 Chiloé event, a moment magnitude of $M_w = 7.48$ is obtained from the inversion of 26 components from 10 seismic stations (9 strong motions and 1 broad-band). The earthquake hypocentre is re-evaluated at 23 km depth with a good RMS of 0.19 and a high confidence of the solution due to the convergence of the possible focal mechanisms.

**Figure S4.** Map of the moment tensor solutions. The vector above each beach ball is the fault geometry (strike, dip, rake). Green circles indicate the centroid of the slip distribution solution of previous studies. Blue circles correspond to the centroid of the moment tensor solutions. Red stars indicate the hypocentre.

**Figure S5.** Nested grid setup used for the modelling of near-field tsunami inundation. White triangles indicate tide gauge stations. (a) Level 1 grid area. Dashed rectangle indicates the level 2 nested grid. (b) Level 2 grid area. Dashed rectangle indicates the level 3 nested grid. (c) Level 3 grid.

**Figure S6.** L-curve of the regularization: The L-curve is a parametric curve, where each point is associated to a lambda parameter. The horizontal axis represents the error between the observed data and the calculated one and the vertical axis represents the solution norm. Higher lambda represents smoother solution, lower lambda
represents smaller error. The optimum regularization parameter is the lambda closest to the corner of the curve. (a) Teleseismic, (b) strong motion, (c) GPS, (d) cGPS, (e) INSAR and (f) tide gauge. The numbers marked are the optimal regularization parameter for each data set, closest to the corner of the curves.

**Figure S7.** Slip distributions for each individual data set inversion. (a) Teleseismic, (b) strong motion, (c) GPS, (d) cGPS, (e) INSAR and (f) tide gauge. Slip scale ranges from blue (0 m) to red (3.8 m). Slip founded for each individual set shows relatively smoother pattern in comparison with the joint inversion.

**Figure S8.** L2 norm of residual versus maximum slip. The red line corresponds to the norm of residual of the slip solution without regularization ($\lambda$ = 0). The goodness of fit is similar to slip distributions with maximum slip between 2.5 and 3.2 m.

**Table S1.** 1-D velocity model used by CSN at the time of the earthquake.

**Table S2.** Magnitude, moment and maximum slip for each inversions.

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.