The Chiloé M_w 7.6 earthquake of 2016 December 25 in Southern Chile and its relation to the M_w 9.5 1960 Valdivia earthquake

Summary

On 2016 December 25, an M_w 7.6 earthquake broke a portion of the Southern Chilean subduction zone south of Chiloé Island, located in the central part of the M_w 9.5 1960 Valdivia earthquake. This region is characterized by repeated earthquakes in 1960 and historical times with very sparse interseismic activity due to the subduction of a young (∼15 Ma), and therefore hot, oceanic plate. We estimate the coseismic slip distribution based on a kinematic finite-fault source model, and through joint inversion of teleseismic body waves and strong motion data. The coseismic slip model yields a total seismic moment of 3.94 × 10²⁰ N·m that occurred over ∼30 s, with the rupture propagating mainly downdip, reaching a peak slip of ∼4.2 m. Regional moment tensor inversion of stronger aftershocks reveals thrust type faulting at depths of the plate interface. The fore- and aftershock seismicity is mostly related to the subduction interface with sparse seismicity in the overriding crust. The 2016 Chiloé event broke a region with increased locking and most likely broke an asperity of the 1960 earthquake. The updip limit of the main event, aftershocks, foreshocks and interseismic activity are spatially similar, located ∼15 km offshore and parallel to Chiloé Islands west coast. The coseismic slip model of the 2016 Chiloé earthquake suggests a peak slip of 4.2 m that locally exceeds the 3.38 m slip deficit that has accumulated since 1960. Therefore, the 2016 Chiloé earthquake possibly released strain that has built up prior to the 1960 Valdivia earthquake.

South America, Earthquake source observations, Seismicity and tectonics, Continental margins: convergent, Subduction zone processes

1 INTRODUCTION

Subduction of the Nazca plate below the South American plate has resulted in various large earthquakes along the Chilean margin (e.g. Lomnitz 2004). The earthquake activity along the Chilean margin is due to the convergence of the Nazca and South American plates, which are colliding at a rate of 6.6 cm yr⁻¹ and an obliquity of 18° in southern Chile (Angermann et al.1999). On 2016 December 25, at 14:22 UTC a magnitude 7.6 subduction earthquake occurred adjacent to the south of Chiloé Island (Fig. 1). In spite of its magnitude and proximity to the city of Quellón, little damage was reported and no tsunami occurred.

Figure 1.

Tectonic setting of the 2016 Chiloé earthquake. The slip distribution of the 1960 earthquake (Moreno et al.2009) is indicated with green lines (5 m slip contours) and plate coupling from Moreno et al. (2011) is shown colour coded. The hypocentre of the 2016 Chiloé earthquake is indicated by a red star. Seismicity is shown colour coded with depth (ISC-GEM catalogue, 1900 January 01–1975 December 31 with stars, gCMT (1976 January 01–2017 February 01 with circles). Two earthquakes from 1919 March are indicated with yellow stars. Oceanic plate ages from Müller et al. (2008). LOFZ after Cembrano et al. (2000). Velocity of the Nazca plate relative to South America after Angermann et al. (1999). Volcanoes are indicated by red triangles.

On the mainland the 1000 km long Liquiñe-Ofqui fault zone (LOFZ, Fig. 1, Cembrano et al.1996) is related to the oblique subduction, and resulting partitioning of deformation along the southern Chilean margin. This partitioning is reflected by the northward movement of a forearc sliver, with sparse current strike-slip faulting along the magmatic arc (Wang et al.2007; Lange et al.2008). The seismicity rate during the last decades (Fig. 2) in Southern Chile is low and seismicity was mostly related to crustal faulting along the LOFZ. There were two phases of seismicity and deformation observed on the LOFZ, which is collocated with the magmatic arc running along the southern Chilean mainland (e.g. Lavenu & Cembrano 1999). The first seismic sequence was related to an M_w 6.2 crustal earthquake in the Aysen region 2007 (45.4°S, e.g. Agurto et al.2012), and was followed by a second sequence related to the eruption of Chaitén volcano in 2008 at 42.5°S (e.g. Watt et al.2009).

Figure 2.

Temporal and trench lateral distribution of seismicity. Top: events from 1920 until 1961 from ISC-GEM Global Instrumental Earthquake Catalogue (Storchak et al.2013) and events from 1961 January until 2017 January 28 from the NEIC catalogue.

The low level of interseismic activity is in contrast with the occurrence of large and great earthquakes related to the plate interface, with a suggested recurrence period of ∼285 yr deduced from the last two millennia (Cisternas et al.2005, 2017). In particular, previous large ruptures occurred in 1575, 1737, 1837 and 1960. The southern Chilean (Valdivia) earthquake of 1960 May 22 is the largest instrumentally recorded earthquake, with a main shock moment magnitude of M_w 9.5 (Cifuentes & Silver 1989). This earthquake initiated at the Arauco peninsula at 38.2°S and progressed ∼1000 km southwards until the rupture terminated near the edge of the subducting Nazca lithosphere, at the intersection of the Chile Rise with the Chilean trench (Fig. 1). This north–south trending rupture produced remarkable changes in land levels (Plafker & Savage 1970). The 1960 main shock was preceded by a series of foreshocks aligned in an NW–SE direction, which started on 1960 May 21 with an M_w 8.1 event (Cifuentes 1989).

The age (and therefore temperature and density) of the subducting Nazca plate decreases drastically from 20 Ma at latitudes of Chiloé Island, to virtually 0 Ma at 46°S where the Chile Rise is currently subducting. The frequency of interseismic activity also decreases substantially toward the Chile Triple Junction (Fig. 2). Estimates for the width of the rupture plane of the 1960 earthquake, which is inclined at an angle of ∼30° below the Chilean mainland (Lange et al.2007; Haberland et al.2009), range between 125 and 150 km (Barrientos & Ward 1990). Thermal models from Völker et al. (2011) show a decreasing width (180 km at 38°S and 80 km at 42°S) of the seismogenic zone (defined here by the 150°C and 350°C isotherms) with increasing latitudes due to the young age (and high temperature) of the subducting oceanic plate.

Plafker & Savage (1970) measured vertical land-level changes at 166 sites along the extent of the 1960 earthquake rupture zone in 1968. Because the vertical displacements were measured eight years after the 1960 earthquake, the data set likely includes post-seismic deformation. Slip models for the 1960 earthquake, inverted from the land-level changes, indicate rupture of the subduction interface along an 850 km long fault, with the majority of slip offshore (Barrientos & Ward 1990). Moreno et al. (2009) inverted the slip of the 1960 earthquake based on a finite-element model using a curved slab geometry, resulting in a smaller amount of slip at larger depths (>80 km), which has previously been interpreted as aseismic slip. Their slip model is characterized by four slip patches with more than 20 m of slip (Fig. 1). One of these slip maxima is located offshore and south of Chiloé Island in the region of the 2016 Chiloé earthquake, midway along the 1960 Valdivia earthquake rupture. In this study, we process seismological data in order to investigate the setting of the 2016 Chiloé earthquake and its relation to the forearc structure, and the 1960 Valdivia earthquake. Based on data from local, regional and teleseismic distances, we locate and determine properties of the main event, fore- and aftershocks of the 2016 Chiloé earthquake.

2 METHODS AND DATA

2.1 Coseismic rupture based on teleseismic observations and local strong motion stations

We used records from eight strong motion stations (three components) located at local and regional distances that were downloaded from the database maintained by the Centro Sismológico Nacional (CSN, www.sismologia.cl, Fig. 3a). The ground-acceleration records were double integrated to obtain displacement, filtered between 0.015 and 0.035 Hz, and resampled to 0.25 s. Broad-band teleseismic waveforms of the FDSN (Federation of Digital Seismograph Networks, Fig. 3b) were retrieved from the IRIS Data Management Centre (www.iris.edu). Teleseismic stations cover an epicentral distance from 30° to 90°. In total 80 P waves (vertical) and 40 SH waves (horizontal transverse) were used. Processing of the teleseismic records included deconvolution of the instrument response, integration to obtain displacement, and windowing 80 s around the body wave arrivals, starting 10 s before the respective P- or SH-wave arrival time. The data were bandpass filtered between 1 and 200 s, and resampled to a 0.25 s sampling rate. The 1-D regional velocity model used by the CSN for Central Chile (Massone, private communication, 2016), was used to compute Green's functions of near-field waveforms and teleseismic body waves. The kinematic finite-source inversion adopted in this study follows the method of Ide & Takeo (1997). Ide & Takeo (1997) expand the slip distribution using 2-D spatial and temporal basis functions, with the expansion coefficients being unknown parameters. Then, the spatiotemporal distribution of slip rate is expanded as a linear combination of basis functions, each one defined by an isosceles triangle in strike, dip and time directions. To stabilize the inversion, temporal and spatial smoothing constraints were used as a priori information. The weighting parameters of smoothing constraints are determined by minimizing the Akaike's Bayesian Information Criterion (Akaike 1980). The non-negative least-squares method of Lawson & Hanson (1974) is used to solve the damped least-squares problem, which ensures positivity of the model parameters. The assumed fault plane strikes at N2°E and follows the slab geometry along dip (Slab 1.0, Hayes et al.2012). This was achieved by subdivision of the fault into four rectangular segments along dip, with widths of 25 km each, dipping with 13°, 15°, 17° and 19°, from top to bottom, respectively. We set a total of 9, 10 and 12 basis functions along strike, dip and time directions, respectively, to represent the slip in one direction. The basis source-time function is an isosceles triangle with duration of 3 s.

Figure 3.

Coseismic slip model of the 2016 M_w 7.6 Chiloé earthquake from joint inversion of teleseismic body waves and strong motion data. (a) Map view showing the distribution of strong motion stations (grey triangles) and the hypocentre of the main shock (red star) located by the CSN. (b) Broad-band stations (blue triangles) at teleseismic distance used in the inversion. Red star indicates the hypocentre of the 2016 Chiloé earthquake. (c) Map view of the final slip distribution. Black arrows correspond to the slip vector scaled to the slip amplitude. Red star indicates the location of the rupture nucleation. (d) Moment-rate function. (e) Fitting of observed (black line) and calculated (red line) ground displacements at the two closest strong motion stations. Station codes and components are indicated on the top of each set of traces. Start time of the records corresponds to the origin time. (f) Comparison of observed (black line) and synthetic (red line) teleseismic body waves at few representative stations. Station code name and the respective P- or SH-wave record is shown inside each box.

We use 12 basis functions in the time domain, with each one having a duration of 3 s and 50 per cent overlap, resulting in the total rise time allowed being 19.5 s. As the focal mechanism is mainly a reverse fault, the slip direction at each gridpoint is allowed to vary between 90° ± 45°. We set the maximum rupture velocity at 92 per cent of the shear wave velocity around the hypocentre depth. This value controls the propagation of a rupture front that expands radially from the hypocentre and defines the onset rupture time of the first knot of the slip rate for each point on the fault when the rupture front reaches it. A slower rupture velocity than the maximum value is allowed using an appropriate choice for the expansion coefficients. In the joint inversion of strong motion and teleseismic data, the relative weighting factors of 1 and 0.02 and were chosen in order to maintain a balance between the error fitting of the teleseismic data in comparison to the measured error of the strong motion data. Note, that these relative weighting factors between the data sets were chosen dependent on the number of available seismograms, effectively upweighting the strong motion stations which are located in the near-field (Fig. 3, panels a and b). The measured error is defined as the L2 norm of the difference between data and synthetics, and normalized by the L2 norm of the data. The relative weighting of data sets is determined prior to including temporal and spatial smoothing constraints. The near-field Green's functions were computed using the numerical code AXITRA (Coutant 1990), which is based on the discrete wavenumber method of Bouchon (1981). Teleseismic body wave Green's functions were computed using the approach of Kikuchi & Kanamori (1991). Synthetic Green's functions were computed with a time step of 0.25 s, and filtered in the same frequency band as the respective observed data sets.

2.2 Fore- and aftershocks

Fore- and aftershocks are based on data from networks C and C1 run by the CSN, and TC and VC networks from Sernageomin (ServicioNacional de Geología y Minería, www.sernageomin.cl). Seismicity was located using NonLinLoc (Lomax et al.2000) and SeisComP3 (www.seiscomp3.org) using a local minimum 1-D velocity model based on amphibious data and optimized for the region of Chiloé Island (Lange et al.2007). Continuous data from up to 54 d prior to the 2016 Chiloé main shock were manually reviewed, and 43 foreshocks based on 872 onset times were located. We reviewed 91 d after the main shock until (2017 March 26) and found 222 locatable aftershocks with 3987 P- and S-arrival times. The mean hypocentral uncertainties are 4.1, 2.7 and 5.6 km for the west–east, north–south directions and depths, respectively (Fig. S1, Supporting Information).

2.3 Regional moment tensor inversion

We inverted regional moment tensors of 17 stronger aftershocks from P and S body waves. For the inversion stations from networks C and C1 from the CSN, TC and VC (Sernageomin), events with epicentral distances of up to 4°, good signal-to-noise ratio (SNR), and local magnitudes larger than M_L 3.8 were included. We used between 6 and 27 stations, and an average of 14 stations. Green's functions were computed using discrete frequency–wavenumber integration (Bouchon 1981), based on a local 1-D velocity (Lange et al.2007), using double-couple point sources.

Before inversion, we deconvolved the instrument response, decimated, rotated and then applied a bandpass Butterworth filter to the observed seismograms in a magnitude-dependent frequency range, predominantly 0.02–0.06 Hz, slightly modified from ranges proposed by other authors (e.g. Kubo et al.2002; Asano et al.2011). Components with low SNR and poor fit were not used for the inversion. Based on the epicentre from the location routine using a local -1D velocity model (Lange et al.2007), we searched for the centroid depth in depth intervals of 2 km. The moment tensor was estimated through a grid search of the best double-couple component in the time domain (Herrmann et al.2011; Herrmann 2013). A grid search over all possible focal mechanisms was performed in increments of 5° in strike, dip and rake. Synthetic seismograms were filtered in the same frequency band as the observed seismograms. First, the waveform fits were evaluated with variance reduction (VR), and then the best solution was chosen based on the largest VR. Fig. S2 in the Supporting Information shows the moment tensor solution and fitting of synthetic and observed waveforms for the largest aftershock.

3 RESULTS

3.1 Coseismic rupture

Fig. 4 shows the final coseismic slip model obtained from joint inversion of teleseismic and strong motion data. The slip distribution is dominated by a large slip patch that covers a region of ∼40 × 30 km². The peak slip reaches ∼4.2 m, located to the north and downdip of the hypocentre, with some slip of about 2 m retrieved to the north and updip of the hypocentre. The total seismic moment computed is 3.94 × 10²⁰ N·m, and gives a moment magnitude of 7.67. The peak moment rate is ∼3 × 10¹⁹ N·m s⁻¹, and occurs 12 s after rupture nucleation, with most of the moment occurring over the first 30 s (Fig. 3d). Figs S3 and S4 in the Supporting Information show the slip models and checkerboard tests when data sets are analysed separately and jointly, respectively. The final slip model obtained from joint inversion of strong motion and teleseismic data explains (1) the deeper slip patch seen when analysing both data sets separately, and (2) the shallow slip (located to the north and updip of the hypocentre) retrieved when only using strong motion data.

Figure 4.

The coseismic slip distribution of the 2016 Chiloé earthquake is shown with a coloured grid and 1 m slip contours. Black dashed lines indicate the depth to the plate interface as defined by the global SLAB1.0 model (Hayes et al.2012). The bathymetry was taken from the GEBCO (https://www.gebco.net/). The red star indicates the hypocentre of the 2016 Chiloé earthquake.

The synthetic seismograms (Fig. 3e) fit well with the observed near-field ground displacements, and the waveform misfit in terms of VR is 89 per cent. We observe that the southernmost strong motion stations (stations A10Y and AY01, Fig. 3a) show a poor fit in comparison to the nearest stations, which we relate to the locally heterogeneous velocity structure. The fit between observed and synthetic teleseismic waveforms is good for both P and SH waves (Fig. 3f). Overall, amplitudes, polarities and the main characteristics of the observed waveforms are well modeled by synthetics. The waveform misfit in terms of VR is 80 per cent. The total VR estimated from the joint inversion of data sets is 84 per cent. Fig. S5 in the Supporting Information shows snapshots of the slip rate taken every 1 s after the rupture nucleation. The rupture propagates rapidly during the first 10 s (reaching the maximum rupture speed), moving mainly to the north and downdip directions. We also observe rupture towards the updip and downdip directions during the first 10 s, but at much less pronounced slip amplitudes. The rupture propagating to the north and in the updip directions slows down after 12 s, and gradually stops. From 12 s onwards, the rupture continues breaking the north and downdip portions of the fault, suggesting a second re-rupturing episode. The overall peak-slip rate retrieved is approximately 1.5 m s⁻¹.

3.2 Hypocentral parameters of seismicity

Interseismic microseismicity from the CSN catalogue (2.1 < M_L < 5.4, 2008 January until 2016 December) and foreshocks (1.3 < M_L < 3.6) occur predominantly in a cluster parallel to the coast of Chiloé Island (Figs 5 and 6), suggesting consistent spatial patterns of seismic moment release. Similarly, during the deployment of an amphibious seismic network in 2005 (Lange et al.2007), seismicity along the plate interface was observed in predominantly thrust type events (M_L 4.0–4.4), with associated aftershock sequences at depths of 12–30 km (white circles, Fig. 5). We observe an epicentral shift of 38 km between the NEIC catalogue (73.951°W, 43.416°S and 35.1 km depth) and our main shock epicentre (74.391°W, 43.517°S, 30 km depth; Fig. 5, yellow symbols). S minus P arrival time difference of the closest station G007 (Fig. 5) is 9.6 s, which is in agreement with a 73 km distance to our hypocentre location. Because NEIC's hypocentre is located at only 40 km distance to the closest station, we think that it is mislocated, presumably due to stronger weighting of stations at regional distances, one-sided station geometry and the use of a global velocity model.

Figure 5.

Fore- and aftershocks of the 2016 Chiloé event together with focal mechanisms from moment tensor inversion. Interseismic activity from 2005 based on an amphibious seismic network (white circles) (Lange et al.2007). The coupling model of Moreno et al. (2011) is indicated with orange lines. The coseismic slip distribution of the 2016 Chiloé earthquake is shown with a coloured grid and 1 m slip contours. High-resolution multibeam bathymetry from R/V Sonne cruise SO181 is encircled with a white line. Hypocentre from NEIC is indicated as a yellow square, hypocentre using the local stations is indicated with a yellow star. Moment tensor for the 2016 Chiloé main shock from gCMT (www.globalcmt.org).

Figure 6.

The 2016 Chiloé earthquake and interseismic activity in the region of Chiloé Island. Circles indicate interseismic activity (CSN catalogue, 2008 January 01 until 2016 December 24 and seismicity from Lange et al.2007) together with the slip model of the 1960 Valdivia earthquake based on the inversion of land-level changes (Moreno et al.2009) shown with green lines (5 m slip contours). Red lines indicate the 1 m slip contours of the 2016 Chiloé earthquake. Volcanoes are indicated with red triangles. Black (CSN) and blue (Sernageomin) triangles indicate local seismic stations used in this study. Blue box indicates location of Fig. 5. The location of the profiles shown in Figs 6 and 7 are shown with a black line and labeled with A–A^΄.

3.3 Moment tensors

Regional moment tensors for aftershocks reveal thrust type mechanisms and centroid depths inline with the geometry of the downgoing slab (Figs 5 and 7). The double-couple mechanisms of the moment tensors indicate faulting at depths of the plate interface. The fault plane of the main shock from the Global Centroid Moment Tensor Catalogue (gCMT) dips 21° with a strike of 7° inline with an easterly dipping plate interface.

Figure 7.

Profiles along the 2016 Chiloé rupture. (a) West–east striking topographic profile. (b) Section showing the aftershock seismicity superimposed on the 2-D vp velocity model from Lange (2008). Red circles indicate aftershocks of the 2016 Chiloé event. Events are shown within 120 km of either side for the profile centred at 43.5°S. The land stations from CSN and active volcanoes are indicated with boxes and triangles, respectively. Model nodes are plotted as crosses. The velocity model is clipped in areas where the resolution of the model is low. See Section 4.4 for explanation of characters ‘W’ to ‘Z’.

4 DISCUSSION

4.1 Coseismic rupture

The total accumulated slip deficit since 1960, based on a plate convergence rate of 66 mm yr⁻¹ (Angermann et al.1999), is 3.66 m for full interplate locking. Taking into account a coupling estimate for the plate interface of ∼90 per cent (Fig. 5 and Moreno et al.2011), the slip deficit is ∼3.3 m. The coseismic peak slip of the 2016 Chiloé event, with 4.2 m, locally exceeds the 3.3 m slip deficit. Undoubtedly, the 2016 Chiloé event locally released slip that has accumulated since 1960. The region where the 2016 coseismic peak slip exceeds the slip deficit is based on two gridpoints in the inversion, and the amplitude of peak slip contains some uncertainties, as indicated by the checkerboard tests (Fig. S4, Supporting Information). If the coseismic slip did indeed exceed the slip deficit that has accumulated since 1960, this implies that the 2016 Chiloé event released slip that accumulated prior to the 1960 Valdivia earthquake (i.e. the 1960 event did not fully release its strain to the same level as the 2016 event). The uncertainties of this slip consideration are inherited from the uncertainty of the 2016 Chiloé coseismic slip model and the coupling model. The effect of slip exceeding the slip deficit becomes even larger when assuming lower coupling values. A coupling of 75 per cent results in a slip deficit of 2.75 m since 1960, of which 1.45 m would have to have accumulated prior to the 1960 Valdivia earthquake. Therefore, the main uncertainties are related to the coseismic slip model and the spatial resolution of the slip and coupling inversions. The peak slip of 4.2 m from our seismological model agrees well with the peak slip inverted from static displacements using GPS receivers, of ∼4.5 m (Ruiz et al.2017) and to the maximum slip of 4.7 m (coseismic and 14 d post-seismic deformation) inverted from InSAR data (Xu 2017). We therefore conclude that the coseismic slip in 2016 may have exceeded the deficit accumulated since the 1960 Valdivia earthquake and the strain release for a given region is larger than expected for slip predictable earthquake recurrence models (which state that the slip of an earthquake is proportional to the time since the prior event (Shimazaki & Nakata 1980). Furthermore, the initial and final stresses of faulting would not be constant over a sequence of earthquakes cycles. In contrast, a peak slip on the order of slip deficit accumulated since 2016 earthquake would be in line with the slip predictable earthquake recurrence models.

The 1960 Valdivia earthquake released an average slip of 17 m along an 850 km long and 130 km wide rupture (Plafker & Savage 1970; Barrientos & Ward 1990), accounting for 250–350 yr of plate convergence (assuming a constant convergence rate). Cisternas et al. (2005) suggest that the 1960 event released the slip that has accumulated since the 1575 event, while two later earthquakes, in 1737 and 1837, left the fault partly loaded with slip deficit that the 1960 earthquake then expended. A similar exceedance of slip compared to the slip deficit since the last large event was described by Nocquet et al. (2017) for the Ecuadorian subduction zone.

The region of the Chiloé 2016 event is characterized by a patch of increased locking (Fig. 1). The almost fully locked zone (>90 per cent) is located slightly updip and south of the coseismic rupture patch (Fig. 5), although the precise location of locking and the shift between peak slip and the locking patch is presumably not well resolved (Moreno et al.2010 and Supporting Information). Moreno et al. (2010) suggest a spatial correlation of interseismic locking with coseismic slip for the M_w 8.8 Maule earthquake in Central Chile. Similar to the Chiloé earthquake, the M_w 8.2 Illapel 2015 earthquake ruptured a region of high locking (Tilmann et al.2016). Our observations from the Chiloé event support the suggestion from Moreno et al. (2010) that interseismic locking might possibly be used to anticipate future ruptures in seismic gaps, given the fundamental assumption that locking and slip are similar. Although the Chiloé 2016 event is spatially located in the rupture zone of the 1960 earthquake, we do not regard this event as an aftershock, due to the 56 yr time difference between these events. The sparse background seismicity (between 41° and 45°S, Fig. 2), and the coupling model (Fig. 1), indicate new strain and stress accumulation during the last decades which were locally released during the 2016 Chiloé event. We cannot completely rule out that the Chiloé 2016 event may have been triggered by relaxation processes within the lower crust and the upper mantle as observed in the 1960 Chile earthquake region south of 38°S (Klotz et al.2001; Moreno et al.2011). For example, using offshore geodetic observations, Sun et al. (2014) observe significant viscoelastic relaxation effects in the aftermath of the 2011 Tohoku earthquake. For the 1960 Valdivia rupture zone westward moving GPS sites south of 38°S, located 300–400 km landward of the 1960 rupture region, suggest a post-seismic mantle stress relaxation or silent slip events on the plate interface at large depths (Khazaradze et al.2002).

4.2 Interseismic activity

Although the background seismicity is very sparse, there are two sequences of historical seismic activity listed in the ISC–GEM catalogue since 1900 (Fig. 2). On 1919 March 2, an M_w 7.2 event occurred at 43°S beneath the South American mainland (Fig. 1, yellow star). This event is listed in the ISC-GEM (Storchaket al.2013) catalogue with a depth of 15 ± 25 km and may therefore be related to the LOFZ. This sequence ended one week later with an M_w 6.8 earthquake further north (Fig. 1). The second phase of pronounced activity is related to outer rise aftershocks (between 43°S and 44.24°S) starting with the M_w 9.5 1960 Valdivia earthquake, and ending in 1965 (Fig. 2 and events indicated by stars, Fig. 1). Later, the background seismicity subsequently returned to a very low rate until the 2016 Chiloé event. Seismicity in the decade before the Chiloé earthquake occurred mostly in a coastal parallel band of events offshore Chiloé Island and parallel to the coastline, likely related to the plate interface (Lange et al.2007).

4.3 Aftershocks

Coseismic slip and aftershocks of the 2016 Chiloé event occurred in the same depth range (Fig. 7), but the 2016 main shock occurred adjacent to, and south of Chiloé Island in a region with very sparse interseismic seismicity. Aftershock magnitudes and locations are not evenly distributed. While larger aftershocks are located in the southwestern aftershock region (Fig. 5), smaller magnitude events are located parallel to the northwest–southeast trending coastline of Chiloé's south coast (Fig. 5). This uneven distribution of aftershock seismicity might be related to the heterogeneous stress distribution along the plate interface. Alternatively, the aftershock distribution may be controlled by deep-reaching NW–SE trending faults which are known elsewhere along the whole South Chilean forearc (e.g. Sernageomin 2003; Glodny et al.2008). However, faults related to the Miocene marine transgressive sediments of Chiloé's south coast are unknown so far (Sernageomin 2003). Furthermore, most of the aftershocks activity is located at depths of the plate interface with sparse seismicity in the overriding crust (Fig. 4).

4.4 Structural setting of the 2016 Chiloé earthquake

The overall crustal structure at latitudes of Chiloé Island (Fig. 7) appears to be very similar to the structure further north at around 38°S (Haberland et al.2009), suggesting limited lateral change along the South Chilean subduction zone. The vp model of the local earthquake tomography study from Lange (2008) reveals features such as the subducting plate (labeled Y, Fig. 7), a wedge-shaped crustal forearc (labeled W, Fig. 7), sedimentary basins in the marine forearc and beneath the longitudinal valley (W and X, Fig. 7). In the western part of the model, the 7.8 km s⁻¹ contour line is inclined following the downgoing slab, whereas in the eastern part of the model the 7.8 km s⁻¹ contour shows a bulge (Z in Fig. 7). A similar, arched 7.8 km s⁻¹vp contour line from local earthquake tomography, was found 500 km further north and interpreted as continental mantle (Haberland et al.2009).

Most of the marine forearc is characterized by very sparse interseismic activity with a sharp onset of events ∼15 km west of Chiloé Island (Fig. 6), at a similar depth to that of the updip limit of the aftershock series of the 2016 event (Fig. 5). This implies that the 2016 earthquake ruptured only the deeper part of the 1960 earthquake rupture zone (Fig. 1). The shallow forearc between the trench and the seismicity band parallel to the coast of Chiloé Island is characterized by an almost complete absence of interseismic seismicity (Fig. 6), but is thought to have ruptured during the 1960 M_w 9.5 Valdivia earthquake (Fig. 1). Here, the plate interface is located beneath the Chiloé basin (Fig. 7b, label W), and the aseismic behaviour of the shallowest few kilometres of the plate interface is inline with global observations of a shallow aseismic zone attributed to underplating of semi-consolidated and unconsolidated sediments (Byrne et al.1988), or to a conditionally stable regime (Bilek et al.2004).

The seismic behaviour in the downdip direction beneath Chiloé Island and the surrounding subduction zone can be classified in the domains suggested by Lay et al. (2012). The Chiloé event ruptured the seismogenic zone (domain B), and left the shallow domain A unruptured (Fig. 8b). Domain A is the region of tsunami earthquakes (Kanamori 1972) that produces large tsunamis relative to their seismic moment (M0), and are usually associated with long rupture duration, for both large and small earthquakes (Bilek et al.2004; Şen et al.2015). The continental crust above the plate interface in the shallow domain A is characterized by a sedimentary basin (Fig. 7 and Scherwath et al.2009), and extends from the trench close to the transition from aseismic to seismic behaviour ∼20 km west of Chiloé Island. A similar seismic band parallel to the coastline, updip of the coseismic rupture and separating the rupture domains in the downdip direction, is observed for some subduction zones such as for the region of the Sumatra M_w 8.7 2005 earthquake (Tilmann et al.2010). Beneath Chiloé Island, both domains A and B broke during the 1960 earthquake, while the 2016 earthquake ruptured only domain B (Fig. 8b). Domain C, located deeper than ∼35 km depth (Fig. 8) did not contain any recorded aftershocks, but is characterized by sparse seismicity during the interseismic period. The location of the 350°C isotherm (based on a shallow inclined slab) at the latitude of Chiloé Island (Oleskevich et al.1999) is close (20 km east) to the intersection of the slab with the continental mantle inferred from the tomography. In general, most of the aftershocks, foreshocks and interseismic activity correlate with temperatures between 150°C and 350°C (Fig. 8a), as expected for the seismogenic zone (Oleskevich et al.1999). For the updip end, Oleskevich et al. (1999) estimate that the plate interface exceeds 100°C 25 km east of the trench. The resulting width of the thermal seismogenic zone at a latitude of 42°S is 100 km. This value is in agreement with the width of the seismogenic zone based on a recent thermal model by Völker et al. (2011), which shows a seismogenic zone based on the 100°C and 350°C isotherms that widens from south to north along the South Chilean margin (from 80 km at 42°S to 180 km at 38°S), due to the older and therefore colder subducting oceanic plate in the north.

Figure 8.

(a) East–west profile showing the distribution of interseismic activity (orange circles: CSN catalogue, 2008 January 01 until 2016 December 24), white circles: interseismic activity from Lange et al. (2007), green circles: foreshock seismicity of the 2016 Chiloé earthquake and red circles: aftershock seismicity of the 2016 Chiloé earthquake. Interseismic activity of the magmatic arc (orange circles ∼72.8°W) is related to the eruption of Chaitén volcano in 2008. Events are shown within 120 km of either side for the profile centred at 43.5°S. The orange lines indicate isotherms from the thermal model of Völker et al.2011. (b) Vertical west–east trending profile for 43.5°S showing the interseismic coupling (Moreno et al.2011), slip of the 1960 Valdivia earthquake (Moreno et al.2009) together with the coseismic slip of the 2016 Chiloé earthquake. Rupture domains following the classification of Lay et al. (2012) are indicated with grey arrows. The land stations and active volcanoes are indicated with boxes and triangles, respectively.

5 CONCLUSIONS

The 2016 Chiloé event is the first significant event within the central and southern segment of the great M_w 9.5 Valdivia earthquake. The 2016 Chiloé event broke a region with locally increased coupling and a large slip of 30 m from the 1960 Valdivia earthquake. The 2016 rupture occurred in a deeper part of the 1960 earthquake, at depths of ∼10–30 km, and in temperature domains inline with the seismogenic zone (150 °C–350°C). The updip limit of coseismic slip and aftershocks of the 2016 Chiloé event is spatially related to interseismic and foreshock activity occurring in a band of seismicity parallel to the coast of Chiloé Island. The faulting style determined using moment tensor inversion of larger aftershocks indicates thrust faulting. Hypocentral depths of foreshock and aftershocks are mostly at depths of the plate interface and the overriding crust is characterized by sparse seismicity. Comparison of the locking model with the coseismic peak slip of the 2016 Chiloé earthquake indicates that strain released in 2016 has accumulated at least since the 1960 Valdivia earthquake. Taking into account the uncertainties of the slip and locking models, we cannot rule out the possibility that some strain released in 2016 might have accumulated prior to the 1960 earthquake. If true, this would imply that the final stress after earthquakes may include stress accumulated over more than one seismic cycle and stress levels reached after the coseismic phase would not be constant over many earthquakes cycles. For the Ecuadorian subduction zone, a similar exceedance of slip in comparison to the slip deficit since the previous earthquake was described by Nocquet et al. (2017). The release of strain that has accumulated at least since 1960 indicates that the 2016 Chiloé earthquake should not be regarded as aftershock of the 1960 event.

Acknowledgements

We thank the Chilean National Seismological Centre (CSN) and the Chilean National Geological y Mining Service (Sernageomin) for providing station data. JR gratefully acknowledges the support from the Chilean National Science Foundation, Project FONDECYT No. 1170804. We are grateful to the Masters and crew of R/V Sonne cruise SO181 and all scientists who were involved in the preparation of the multibeam bathymetry. We thank Jess Hillman who contributed to this manuscript by carefully spell checking and correcting grammar. Our special thanks go to Gavin Hayes and an anonymous reviewer for their constructive comments and suggestions.

SUPPORTING INFORMATION

Supplementary data are available at GJI online.

Figure S1. Mean uncertainties based on error ellipsoids (68 per cent confidence) from NonLinLocfor 0.05° × 0.0.5° tiles for the fore- and aftershocks. (a) mean longitude uncertainty, (b) latitude uncertainty and (c) depth uncertainty.

Figure S2. Moment tensor inversion of the largest aftershock of the Chiloé 2016 event (M_w 5.5, 2017 January 03 21:19:08 UTC). (a) Synthetic (blue lines) and observed (red lines) waveforms for the best centroid depth for selected stations used for the inversion. Each component is labeled with its maximum amplitude (up left), time-shift (in seconds) and data fit in percent (up right). (b) Waveform fit as a function of source depth showing best moment tensor solution for each depth. (c) Best double-couple solution. Please note, that Regional moment tensor solution for this aftershock is quite similar to gCMT catalogue (there is currently (2016 January 28) no other gCMT solution for the aftershock sequence available). The parameters from the gCMT catalogue are similar (strike=163°, dip=70°, slip=93°, M_w = 5.5, M_s = 5.5 and M₀ = 2.14e + 24 N·m).

Figure S3. Coseismic slip models obtained from kinematic finite-source inversion using (a) teleseismic body waves, (b) strong motion and (c) both teleseismic and strong motion data sets. (d) Comparison of moment-rate functions obtained when analysing each data set separately and jointly. Grey stars in each slip model represent the hypocentre location. Black arrows correspond to the slip direction scaled to the slip amplitude. The colour bar on the left of each fault model indicates the slip amplitude in metres.

Figure S4. Coseismic slip models obtained from the resolution test. (a) Input slip model. (b) Slip obtained from the inversion of teleseismic data. (c) Slip retrieved when using only strong motion data. (d) Slip distribution computed from joint inversion of both the teleseismic body waves and near-field data. Grey star indicates the hypocentre. The arrows indicate the slip vector scaled to the slip amplitude. We used the same fault model, stations and parameter configuration used in the inversion of the coseismic slip model of the main shock. The input synthetic slip model shown in the figure (top left) was used to generate synthetics waveforms at teleseismic and strong motion stations. We added random Gaussian noise having a standard deviation of 5 per cent of the maximum amplitude to each synthetic waveform. The local slip-rate time function using in the forward modeling corresponds to a single isosceles triangle with a duration of 15 s. The input faulting mechanism is a pure reverse fault, but during the inversion the slip is allowed to vary between 45° and 135°. Following the inversion method described in our study, we run three inversions using (1) only teleseismic data (panel top right), (2) only strong motion (bottom left) and (3) both, teleseismic and strong motion data(bottom right). The results show that the slip pattern retrieved using joint inversion performs better than using each data set separately.

Figure S5. Snapshots of the slip rate taken every 1 s after the rupture initiation (elapsed time is indicated at the top of each fault plane). Yellow contours correspond to the rupture front that propagates at the faster rupture speed set in the inversion. Colour bar indicates the slip rate amplitude in m s⁻¹.

Please note: Oxford University Press are 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.