Modelling cascading ruptures on near-orthogonal strike-slip fault system: the 2019 Cotabato (the Philippines) earthquake sequence

Earthquake sequences on near orthogonal strike-slip faults are not uncommon, as observed both in the subduction zone outer-rise region (e.g., the 2000 and 2012 Wharton basin earthquakes off-Sumatra, Robinson et al., 2001; Wei et al., 2013) and in the shallow continental crust (e.g., the 1987 Superstition Hills earthquakes, Hudnut et al., 1989; the 2019 Ridgecrest earthquakes, Shi and Wei, 2020). However, according to the Coulomb faulting theory (Anderson, 1951), strike-slip faults intersecting at an angle of around 90° (at 45° from the maximum principal stress) would require a near-zero friction coefficient, which is not consistent with the observed values 0.6~1 in nature, e.g., deep borehole stress measurements (Townend, 2006). Thus, the mechanisms controlling these cascading orthogonal ruptures remain poorly understood. The 2019 Cotabato earthquakes (the Philippines) provide a new opportunity to further explore the driving mechanism causing orthogonal strike-slip earthquakes, since abundant geodetic and seismic data sets are recorded in this sequence.

In this research, we focus on four Mw6.4+ events during the 2019 Cotabato earthquake sequence. Since the fault geometry is critical to analyze the potential stress triggering between earthquakes, careful processing and modelling of the data sets are required to provide a robust and reliable fault geometry. To better constrain the fault geometry, two types of observations were utilized, surface displacement data with a high spatial resolution and ground motion data with a high temporal resolution. (a) Geodetic modelling. We acquired eight ALOS-2 L-band SLC images, and generated ten interferograms monitoring the ground displacement, including seven ascending and three descending interferograms. To avoid the influence of phase unwrapping errors, we improved and applied an art-of-the-state Bayesian geodetic inversion approach (Jiang and González, 2020) by using the InSAR wrapped phase and allowing the estimation of multiple fault geometry simultaneously. (b) Seismic modelling. Seismic waveforms were collected from IRIS, including one regional station and over ten teleseismic stations. We performed Multiple Point Source inversions (Shi et al., 2018) to determine the subevents’ location and double-couple focal mechanism. Geodetic and seismic inversion results were cross-verified and updated to reconcile both datasets. Our results show that (1) in mid-October, the first Mw6.4 earthquake occurred on an NW-SE-striking fault at the depth range of 10-18 km; (2) the second Mw6.6 earthquake ruptured the shallow part of the same fault, followed by the third Mw6.5 earthquake two days later but rupturing a NE-SW-striking fault; (3) in mid-December, the most energetic Mw6.8 earthquake occurred on an NW-SE-striking fault, located at SE of the first two events. Coulomb stress analysis suggests that the friction coefficient on the NE-SW-striking fault has to be very low to allow the rupture on the near orthogonal faults. Our results indicate that the earthquake sequence is a cascading rupture that involved both weak and strong faults in which pore fluid pressure may have played a key role.