Did the September 2010 (Darfield) earthquake trigger the February 2011 (Christchurch) event?

We have investigated the possible cause-and-effect relationship due to stress transfer between two earthquakes that occurred near Christchurch, New Zealand, in September 2010 and in February 2011. The Mw 7.1 Darfield (Canterbury) event took place along a previously unrecognized fault. The Mw 6.3 Christchurch earthquake, generated by a thrust fault, occurred approximately five months later, 6 km south-east of Christchurch's city center. We have first measured the surface displacement field to retrieve the geometries of the two seismic sources and the slip distribution. In order to assess whether the first earthquake increased the likelihood of occurrence of a second earthquake, we compute the Coulomb Failure Function (CFF). We find that the maximum CFF increase over the second fault plane is reached exactly around the hypocenter of the second earthquake. In this respect, we may conclude that the Darfield earthquake contributed to promote the rupture of the Christchurch fault.

Earthquakes interact. The research in the last decades demonstrated that over major active faults or fault systems, where seismologists registered the occurrence of an earthquake, the probability of occurrence of a second shock increases or decreases according to stress changes1234. Indeed, a mainshock perturbs the stress state in other sections of the same fault or in adjacent faults: this theory is known as Coulomb Stress Triggering5. The hypothesis is that once an earthquake occurs, the stress does not dissipate, but it propagates in the surrounding area, where it may increase the probability of occurrence of further earthquakes. Several examples can be found in the literature: in some cases, two seismic events, Earthquake-1 and Earthquake-2, occur over a long time span4; in other cases they are temporally very close678. The towns of Darfield (Canterbury) and Christchurch have been hit by strong earthquakes within a time span of a few months9101112. Earthquake-1 took place on September 03, 2010 (at 16:35:46 UTC, Mw 7.1), while Earthquake-2 occurred on February 21, 2011 (at 23:51:43 UTC, Mw 6.3), a few kilometers more to the East. Between the two events, more than 4.000 aftershocks were recorded over the Greendale fault and other secondary structures. Could the second event have been triggered by a stress redistribution induced by the first one? The answer is in terms of probability change. Indeed stress change can trigger a second event if fault conditions are close to failure, otherwise it will hasten the occurrence of future earthquakes.

New Zealand is located across the margin of the Australian and Pacific plates in the southern pacific. Here the relative obliquely convergent plate motion varies from 30 mm/yr in the south of the country to about 50 mm/yr in the northern part13. Starting from the Middle-Late Cretaceous and the Late Cretaceous-Paleocene episodes of rifting, New Zealand separated from Australia and Antarctica14 since the Late Paleocene, during which a new extensional plate boundary initiated about 45 Ma15. Around 25 Ma the processes of oblique compression across the plate boundary initiated16, with an increase of convergence rates since the late Miocene.

The study area is located within the Canterbury-Chathamas platform, characterized by a 1000-km-long dextral transpressive zone, represented by the Alpine fault and the Malborough fault belt (Figure 1a) 17. This fault accommodates the relative plate motion between the NW-dipping Hikurangi and SE-dipping Puysegur subduction zones18. On the Alpine fault, representative strike-slip rates are 25–30 mm/yr1920 and only geological data provide evidence of large-to-great earthquake occurrence, with recurrence times of hundreds of years. The slip rates observed on the Malborough fault belt are greater than 1 mm/yr; here many moderate-to-large earthquakes occurred in historic times, including the Mw 8.1–8.2 1855 Wairarapa and the Mw 7.8 1931 Hawkes Bay earthquakes. This strike-slip province is characterized by several fault strands, hundreds of kilometers in length, associated with folds and reverse faults. Very few data are available about the distribution of segment boundaries on the main strands and to constrain the timing of paleoseismic events.

Figure 1

(A) Tectonic setting of southern New Zealand island 14.

Black lines indicate the main tectonic structures discussed in the text. Red rectangle represent the study area. (B) Aftershocks distribution between the two main earthquakes on the September 3, 2010 and on February 21, 2011. Red lines represent the activated faults. (C) Unwrapped interferograms in LOS geometry. Red stars are the hypocenter location (from CMT) of the two events; Focal mechanisms are also shown.

The Mw 7.1 Darfield (Canterbury) earthquake occurred along a previously unrecognized east-west fault line, the strike-slip Greendale fault (Figure 1b). In this area no active faults had previously been mapped, nor are large historical earthquakes known to have occurred. The only active faults known in this region are located further west, at the foothills of the Southern Alps, where several M > 6–7 earthquakes have occurred in the past 150 years. The teleseismic moment tensor (USGS) and finite fault solutions2122, providing a far-field observation of the earthquake, have indicated a dextral strike-slip fault (Figure 1c) in agreement with both the orientation and the sense of slip of the documented surface rupture11. On the other hand, the near field seismological observations (first motion and regional moment tensor,10) show a large reverse faulting component, in contrast with the teleseismic solutions. This is due to differences in the two measurement techniques, which analyze high (near field) and low (teleseismic) frequency waves and are therefore sensitive to small or large features respectively. A possible interpretation is that the earthquake started as a thrust event in a smaller scale structure and continued over the main structure (Greendale fault) with a strike-slip mechanism, as indicated by the teleseismic solutions, thus accommodating the regional stress and releasing the largest energy fraction.

The 6.3 Mw Christchurch earthquake occurred on February 21st, 2011, approximately 5 months after the Darfield earthquake. The hypocentre was approximately 6 km south-east of Christchurch's city center, at a depth of 5–6 km, generated by a blind ENE oblique-thrust, faulting at the easternmost limit of the Darfield aftershocks (Figure 1b) 12. No specific structure for this event is directly linked to the main fault of the 2010 mainshock. The focal mechanism solution indicated a right-lateral oblique thrust faulting mechanism (Figure 1c).

Roughly 20 years ago the Earth Sciences received the impact of the newborn Differential Synthetic Aperture Radar Interferometry (DInSAR) technique23. DInSAR has become a key element in multidisciplinary studies of earthquake24. In order to investigate Earthquake-1 and Earthquake-2 we apply DInSAR, using two pairs of satellite images acquired by the Japanese mission ALOS (Advanced Land Observing Satellite), and its onboard PALSAR (Phased Array type L-band Synthetic Aperture Radar) sensor, along two adjacent tracks. We use the Shuttle Radar Topographic Mission digital elevation model25 to remove topographic fringes from the interferograms (Figure 1c).

Results

The map of the surface coseismic displacement (Supplementary Figure S1) of Earthquake-1 shows a complex pattern, meaning that the Greendale fault and other secondary buried faults moved during the seismic event. Although this is not an evidence of the triggering of Earthquake-2, it can be noted that the surface displacement field extends near to the epicenter of the February 21st earthquake. The latter occurred about 50 km East in a less complex scenario; we argue, based on DInSAR, that only a single fault is rupturing during the shock12 (Supplementary Figure S1). Looking at both interferograms it is plain that PALSAR is an effective tool to capture the surface displacement field in case of moderate and strong earthquakes2627. We apply an adaptive filter28 in order to reduce the noise, and a Minimum Cost Flow (MCF) phase unwrapping algorithm29.

Based on preliminary field observations11, the Greendale rupture zone is unequivocally the area of main deformation, with significant amounts of displacement. The average displacement is ∼2.5m along the ∼30km long main rupture and reaches a maximum of ∼5m suggesting that the greatest energy fraction is released by the strike-slip Greendale fault. Smaller scale ruptures (blind thrusts) on nearby faults are also associated with this earthquake, but are expected to contribute less to the CFF estimation.

The slip distribution for the Darfield earthquake is estimated by linearly inverting the sub sampled DInSAR deformation map. In a previous work9 the authors have investigated the coseismic deformation produced by the Darfield earthquake, following the main surface rupture, in order to constrain the Greendale fault and its related extensions. We fix our main fault in a similar way, following the documented surface rupture. We approximate the Greendale fault as a concatenation of three planar rectangular strike-slip fault planes, with 44km of total length and 12km of width. These are the main segment striking E-W and coincident with the major part of the mapped surface rupture, the NE-SW extension of the fault and the step-over (offset to the north) segment at the eastern end of the Greendale fault. The fault planes are subdivided into a discrete number of rectangular patches30. A Green's function matrix (kernel) is composed by imposing a unitary dislocation on each patch and subsequently collecting the E, N and vertical surface displacement components projected along the satellites line of sight (LOS). The linear inversion is performed adopting the Occam's smoothing scheme31, minimizing the chi-square and the second order derivative (Laplacian) to avoid large, unphysical oscillations in slip values. During the inversion process we also take into account three thrusts: the first one is near the hypocenter, the second at the NW end of the Greendale fault, distant from the main deformation field, and the third near the step-over segment. The largest part of the slip is concentrated over the three segments of the main fault, with a peak value of 6.5m in the central one, suggesting that the majority of the energy is released by the Greendale fault (Figure 2). The resulting geodetic moment is ∼5×1019 Nm (Mw ∼7.1) if a crustal rigidity of 30GPa is assumed in accordance with seismological estimates32. Our results for the slip distribution reveal a similar pattern to the one derived in9.

Figure 2

(Top view) Selected frame of the East coast of South New Zealand.

Surface projection of the fault planes adopted for the Darfield 2010 (Earthquake-1) and for the Christchurch 2011 (Earthquake-2) events. The sense of the slip for Earhquake-1 is right-lateral (180°) for the Greendale fault and its lateral extensions, as indicated by the red arrows, and thrust (90°) for the three minor faults. The black arrows indicate the dip direction of the three minor fault planes. Aftershocks distribution (black dots) corresponds to the time period between September 3rd, 2010 and February 23th, 2011 12. (Bottom panel) A 3D perspective view for Earthquake-1 and Earthquake 2. The largest part of the slip (max ∼6.5m) is concentrated in the middle segment (Greendale fault) from 0 to 6km depth. Coulomb stress change is estimated for the Earthquake-2 fault plane. The red and black stars indicate the hypocenter of Earthquake-1 and Earthquake-2 respectively (GNS Science). Both panels are in UTM WGS84 coordinate system.

In order to investigate the consistency of our results, we performed an uncertainty analysis, computing the covariance and resolution matrices for our model parameters33. The analysis was limited to the three main segments of the Greendale fault, for the sake of simplicity. The resolution matrix (R) can be used to estimate the spatial resolution at the location of any single fault patch in the model. If a single slip patch is perfectly resolved, the corresponding column in R will have a value of 1 at the main diagonal position, and the off diagonal elements will be zero. Results indicate (Supplementary Figure 3) that only the first row of patches is well resolved (R∼0.9) and that the resolution decreases rapidly below 2km depth. The results for resolution indicate that only general slip features can be resolved at depth. The standard deviation is obtained from the square root of the elements on the main diagonal of the covariance matrix. The uncertainty bounds are acceptable and range between 0 and 70cm everywhere in the fault plane, with the exception of the Western fault segment, where the error reaches higher values at depth (SE corner).

The DInSAR technique offers a useful tool for fault characterization34. Indeed the fringe shape, rate and orientation can be related to fault parameters, such as geometric dimensions and orientation angles. We have used a novel approach35, based on the Okada model30 and Neural Networks (NNs), to investigate the fault geometry of the Christchurch earthquake. One of the advantages of this method is that it rapidly achieves a determination of the rupture plane.

Once the geometries of the two faults are defined (see Table 1), we focus our analysis to understand the role of the first earthquake in promoting the rupture of the second event through the evaluation of the Coulomb Failure Function (CFF). The CFF is obtained by computing the stress tensor corresponding to the elastic dislocation induced by the Canterbury earthquake, projecting it onto the rupture plane of the Christchurch earthquake, and evaluating the relative weights of normal and shear stresses, assuming an effective friction coefficient of 0.4. This value is in agreement with laboratory values of friction and moderate pore pressure when fluids are not fully expelled4536. Positive or negative variations of the CFF indicate that the stress field is acting to promote, or oppose, the rupture, respectively. Our results (Fig. 3a) show that the rupture of the 2010 Darfield earthquake loaded a large portion of the crust with stress values exceeding 1 bar. If we take into account the three-dimensional location of the Christchurch rupture plane (Fig. 2), we see that that the shallower part of the fault (down to about 5km depth) has actually been unloaded by the Canterbury earthquake. On the contrary, the remaining, deeper part of the fault has been brought closer to rupture, with largest values of stress loading towards the south-western tip of the plane. The average CFF value on the loaded portion of the fault is over 0.01 MPa, with peak values exceeding 0.03 MPa. Remarkably, the peak stress increase occurs in the southwestern part of the fault, where rupture nucleation as occurred according to GNS localizations. These stress levels are definitely non-negligible, since a stress value of the order of 0.01 MPa is considered a threshold for effective triggering of seismic events37. We expect these estimates to be affected by several uncertainties. First, the geometric parameters of the reconstructed planes are known within a certain level of precision, depending on the DinSAR data coverage and technical details of the modeling procedure. Also the assumption of an elastic half-space is likely to introduce a bias, given the tectonic framework of the region. Finally, postseismic deformations add a time-dependent stress perturbation that contributes to the total effect on the fault. Quantifying the impact of all these contributions turns out to be a quite complex task; however, they are not expected to alter the qualitative aspect of our conclusions. In Fig. 3b we compare the location of aftershocks with CFF variations on optimally oriented fault planes, assuming a background regional stress corresponding to an E-W compressional tectonics3839. Even if a complete analysis of the sequence seismicity would be beyond the scope of this work, we see from Fig. 3b that most of the aftershocks occur in loaded regions, suggesting the hypothesis that the sequence evolution may indeed be driven by stress redistribution mechanisms.

Table 1

The parameter retrieval of the Earthquake-1 and Earthquake-2, obtained by means of Okada model and NNs.

Earthquake-1
Fault Segments	Length (km)	Width (km)	Depth (km)	Dip angle (deg)	Strike angle (deg)	Rake angle (deg)	Mean Slip (cm)	Lat (deg) North (km)	Lon (deg) East (km)
1st segment NE-SW	12	12	0	87	269	180		−43.5817 (5173600.0)	172.4184 (614520.00)
2nd segment Greendale Fault	20	12	0	87	265	180		−43.5886 (5173000.0)	172.2969 (604700.00)
3rd segment Step Over	12	12	0	87	302	180		−43.599 (5172040.0)	172.0655 (586000.00)
4th segment (thrust near step over)	8	8	1	65	24	90		−43.576 (5174400.0)	172.278 (603220.00)
5th segment (thrust near hypocenter)	8	8	1	75	31	90		−43.5938 (5172650.0)	172.100 (588800.00)
6th segment (thrust near NW-SE extension)	8	8	0.5	65	214	90		−43.5471 (5177979.0)	171.951 (576891.69)

Figure 3

(a) Coulomb stress change projected onto the Earthquake-2 rupture geometry estimated on a horizontal plane at 6 km depth.

A black rectangle marks the surface projection of the 2011 Christchurch earthquake (Earthquake-2). Yellow stars mark the location of the two events (GNS Science). (b) Coulomb stress change induced by Earthquake-1 on optimally oriented fault planes. Circles mark the epicentral location of aftershocks from the GNS catalogue.

Discussion

The DInSAR results allowed us to investigate a relation between a pair of spatio-temporally close earthquakes. The interferograms provided a first input to our analysis chain, in the form of surface displacement fields to be used for inversion modeling. Subsequently, using the fault geometries of both earthquakes and the slip distribution of the first one (Darfield earthquake, September 3, 2010), we calculated the CFF over the second fault plane.

After the September 3 earthquake, a sequence began with aftershocks located along the Greendale fault and some hidden secondary faults. A magnitude 5.1 aftershock occurred on September 8 nearby the epicenter of the February 21 event. Even though this event occurred five months after the September 3 earthquake, some scientists consider the February 21 an aftershock caused by a fault rupture within the zone of aftershocks of the mainshock40. Research on earthquake triggering investigated earthquake interaction at different scales: mainshock-mainshock and mainshock-aftershock (41 and references therein). Furthermore, it is now accepted that earthquake triggering occurs at all scales and that there is no mechanistic difference between the origin of foreshocks, mainshocks and aftershocks. Based on such premises, Earthquake-2 can be interpreted as a mainshock along a second fault, promoted by the stress perturbation of Earthquake-1, or as the largest aftershock of the sequence started with Earthquake-1. It is however considered out of the scope of the present work to pursue this issue further, although it would be interesting from the statistical point of view. With our analysis we cannot state whether Earthquake-2 was (or was not) extremely unlikely to occur40 without the preceding Earthquake-1; however, the outcome of our work is that Earthquake-1 has loaded the Earthquake-2 fault, bringing it closer to failure.

Methods

We selected a dataset of images from the Japanese ALOS satellite. ALOS has onboard a Synthetic Aperture Radar sensor, PALSAR, which is an active microwave sensor using L-band frequency to achieve cloud-free and day-and-night land observation. ALOS PALSAR allows a Fine Beam Single (FBS) polarization mode that achieves a spatial resolution of 5 m in azimuth and 7 to 44 m in ground range, according to the incidence angle, and a Fine Beam Dual (FBD) polarization mode, which, compared to the FBS, has the same azimuth resolution and a halved range resolution (14 to 88 m).

In this study we use two pairs of ALOS PALSAR images along two adjacent ground tracks. The image pair relative to the September 3rd, 2010, earthquake is in FBD mode with an incidence angle of about 38°. The pre-seismic image dates August 13, 2010 and the post-seismic September 25, 2010. The image pair used to study the second earthquake, February 21, 2011, is in FBS mode with an incidence angle of about 38°. In this case pre- and post-seismic data are taken on January 1, 2011 and February 25, 2011, respectively.

We apply the DInSAR technique to detect coseismic ground deformation. The interferograms have been computed with a square pixel of about 28 m achieved by applying a 3 by 8 multi-look factor in slant-range and azimuth respectively for FBD, and a 6 by 8 factor for FBS, in order to improve the signal-to-noise ratio. Furthermore, the 90-m Shuttle Radar Topography Mission (SRTM) digital elevation model25 has been used to simulate and remove the topographic contribution of the interferometric phase. Both interferograms maintain a good coherence, which allows capturing most of the coseismic deformation, also for the second earthquake where a lot of damages occurred in the city, confirming the effectiveness of the PALSAR L-band sensor in the case of strong events27. In order to mitigate phase noise we apply an adaptive filter28, while to retrieve the Line Of Sight (LOS) displacements, we use a minimum cost flow phase unwrapping algorithm29.

After mapping the coseismic ground deformation we perform a linear inversion in order to retrieve the slip distribution over the Greendale fault. We take into account the documented surface rupture by representing the Greendale fault as a concatenation of three smaller faults. The geometric features of the rectangular planes are as follow: -1- central segment (Greendale fault), E-W orientation, Length 20km, Width 12km, dip 87° -2- western segment (NW-SE extension), NW-SE orientation, Length 12km, Width 12km, dip 87° -3- eastern segment (step over toward north), E-W orientation, Length 12km, Width 12km, dip 87°, -4- thrust fault near step over eastern segment, Length 8km, Width 8km, dip 65°, -5- thrust fault near hypocenter, Length 8km, Width 8km, dip 75°, -6- thrust fault near the NW-SE Greendale fault extension, Length 8km, Width 8km, dip 65°. The rake vector is fixed at 180°(right lateral faulting) and 90°(thrust faulting) for -1-, -2-, -3- and -4-, -5-, -6- respectively (Supplementary Figure S4). Subsequently the fault model is subdivided into square patches of 2km per side. We compose the Green's function matrix (kernel) by imposing a unitary dislocation on each one of the patches, calculating the E, N and vertical surface displacement components, and projecting along the satellites LOS.

In order to retrieve the fault geometry of the Christchurch earthquake we use an innovative approach based on Neural Networks35, that calculates the surface displacement field due to a dislocation on the fault plane at depth assuming an elastic half space, and using the analytical solutions provided by Okada30. This method has been chosen to investigate this second earthquake because it provides the major advantage of allowing a fast development of a preliminary model for the seismic source and thus a fast determination of the rupture plane for the Christchurch event. Moreover, it avoids the requirement for phase unwrapping.

Neural networks (NNs) have already been recognized as being a powerful tool for inversion procedure in remote sensing applications. They are composed of an ensemble of nonlinear computational elements (called neurons) connected by the so-called synapses, each characterized by a synaptic weight. Compared to other methods, the use of NNs is often effective, because they can simultaneously address nonlinear dependences and complex physical behavior with reduced computational efforts and without the need of any a priori information.

We then investigate the role of the first earthquake in promoting the rupture of the second event through the evaluation of the Coulomb Failure Function (CFF). First we compute the elastic strain tensor corresponding to the dislocation field induced by the Darfield (Canterbury) earthquake, using the analytical solutions provided by Okada30. By applying the relations of standard elasticity theory, the spatially-variable strain is converted into an incremental stress tensor that acts as a perturbation to the pre-existing (unknown) stress state of the crust. Once the perturbation to the stress tensor is known, its effect on a given fault mechanism can be assessed by evaluating the CFF variation, defined as ΔCFF = Δτ + μ(Δσn + ΔP), where Δτ and Δσn are, respectively, the shear and normal stress changes, μ is the friction coefficient and ΔP is the pore pressure change37. It is common to rewrite the definition of the CFF variation as ΔCFF = Δτ +µ'Δσn, where µ' is an effective friction coefficient that takes into account static friction, hydrostatic pressure and pore fluid pressure45; in our work we assumed µ' = 0.4. The knowledge of ΔCFF allows to establish whether the imposed stress field acts to promote (ΔCFF>0) or to oppose (ΔCFF<0) the dislocation on the considered fault. As a rule of thumb, in literature a CFF increase of 0.01 MPa (the magnitude of tidal loading) is considered a threshold for the effectiveness of triggering 45.

Author Contributions

S.S., M.C. and C.B. contributed with the DInSAR processing. C.K, M.M., M.S. provided the geological and tectonic analysis. C.K. and M.P. have performed the modeling inversion. D.M. provided the CFF. All authors contributed to write the manuscript and to discuss the results.