Shallow Megathrust Rupture during the 10 February 2021 $Mw$ 7.7 Southeast Loyalty Islands Earthquake Sequence

The Vanuatu (New Hebrides) subduction zone extending from about 10.5° S near the Santa Cruz Islands to about 23° S southeast of the Loyalty Islands is a seismically active arc with the Australian plate underthrusting the Pacific plate and the western and southern North Fiji basin microplates (e.g., Louat and Pelletier, 1989; Fig. 1). The arc is laterally bounded by strong corner bends, with transforms connecting to the Solomon Islands in the north and to a seismic line south of Fiji in the south. The northern Vanuatu subduction zone has experienced larger historical earthquakes than the southernmost portion (McCann et al., 1979; Nishenko, 1991; Cleveland et al., 2014), with the most recent being the great $Mw$ 8.0 event of 6 February 2013 near the Santa Cruz Islands (Lay et al., 2013; Cleveland et al., 2014; Hayes et al., 2014). Large earthquake doublets and triplets have occurred both in northern Vanuatu and the Loyalty Islands region, indicating strong stress change interactions along the megathrust plate boundary (e.g., Cleveland et al., 2014).

The southern bend in the subduction zone near the Loyalty Islands intrinsically involves complex tearing and deformation of the underthrusting Australian plate, with arc‐normal thrusting occurring all along the ∼70° bend in the zone. This curved geometry is accommodated by tearing the underthrusting plate, with an underthrust flange along the eastern extension sustaining convergence as there is a transition to strike‐slip faulting to the east (e.g., Monzier et al., 1984; Lister et al., 2012). Strike‐slip shearing deformation within the overriding southern North Fiji basin block appears to help accommodate the progressive slip partitioning around the bend (Louat and Pelletier, 1989). Observations from the Global Positioning System along the Vanuatu archipelago and the Loyalty Islands (e.g., Calmant et al., 2003) show variable convergence along the western North Fiji basin complicated by deformation and segmentation into blocks within the overriding plate. In the southern region, the incoming plate has significant bathymetric relief, with the Loyalty ridge impinging on the subduction zone.

This region has a substantial history of major earthquakes (Fig. 1a), with magnitude as large as $Mw$ 7.7, but no history of great (⁠$Mw≥8$⁠) earthquakes. The convergence rate is high, ∼8 cm/yr, and the background seismicity is correspondingly high, but the region is distinct in having an unusual amount of activity seaward of the trench. Large earthquakes in this region are accompanied by highly productive aftershock sequences. Numerous sequences of large historic ruptures along the southern Vanuatu arc and the Loyalty Islands (Fig. 1a) show strong interactions between the megathrust and the intraplate environment seaward of the trench. The 1995 $Mw$ 7.7, 2004 $Mw$ 7.1, and 2018 $Mw$ 7.5 events west of the trench all involve normal faulting, whereas the 1980 $Mw$ 7.5, 2003 $Mw$ 7.3, and 2018 $Mw$ 7.1 events all involve thrust faulting. The 1980 $Mw$ 7.5 event was part of a sequence of $Mw$ 6.5–7.5 thrust events (Vidale and Kanamori, 1983), and the 1995 $Mw$ 7.7 event was an isolated normal faulting event relatively far seaward of the trench, whereas the 2003–2004 and 2018 events are doublets, closely spaced in time and location.

On 10 February 2021, the plate boundary southeast of the Loyalty Islands experienced a major thrusting earthquake with $Mw$ 7.7 (13:19:55.53 UTC, 23.051° S, 171.657° E, 10.0 km deep; U.S. Geological Survey National Earthquake Information Center [USGS‐NEIC]) (Fig. 1a). The earthquake sequence was very productive, including three large thrusting foreshocks with $Mw$ from 5.6 to 6.1 within 1 hr before the mainshock and 347 aftershocks larger than $Mw$ 4.5 in the first 45 days after the event (Fig. 1c; Ye et al., 2020; Fan et al., 2021). We analyze the rupture process of the mainshock using global broadband seismic recordings and consider the 2021 sequence relative to other sequences that have occurred along the Loyalty Islands.

The 10 February 2021 $Mw$ 7.7 southeast Loyalty Islands mainshock involved thrust faulting on a northward dipping fault with a strike close to east–west, centered in the strongly curved plate boundary between the Australian plate and southern North Fiji basin. The USGS‐NEIC W‐phase inversion has a predominantly double‐couple solution with strike $ϕ=246°$⁠, dip $δ=17°$⁠, rake $λ=65°$⁠, and a seismic moment $M0=4.364×1020N·m$ (⁠$Mww$ 7.69). The centroid depth is 25.5 km. The Global Centroid Moment Tensor (Global CMT) solution, based on a dataset of 467 body waves, 517 mantle waves, and 440 surface waves low‐pass filtered at periods of 50, 150, and 50 s, respectively, has a best double couple with $ϕ=277°$⁠, $δ=22°$⁠, $λ=98°$⁠, $M0=4.12×1020N·m$ (⁠$Mw$ 7.7), centroid depth 21.0 km, and 28.5 s centroid time shift. Given the relatively large differences in strike and rake between these routine solutions, we performed our own W‐phase inversion (Kanamori and Rivera, 2008) using 145 ground‐motion channels from 79 stations filtered in the passband 0.002–0.005 Hz. We carefully screened the data quality and checked the waveform predictions. This inversion gives a centroid depth of 25.5 km, the best double couple having $ϕ=276.9°$⁠, $δ=21.2°$⁠, $λ=94.8°$⁠, $M0=4.15×1020N·m$ (⁠$Mw$ 7.68), and a 27 s centroid time shift, which is very similar to the Global CMT solution.

The USGS computed an initial planar finite‐fault solution for the $Mw$ 7.7 southeast Loyalty Islands event using teleseismic body and surface waves, assuming a strike of 264° and a dip of 19°. The result is a concentrated slip patch at depths less than 10 km with slip of up to 6 m. In our own preliminary inversions of teleseismic body waves following the procedure of Ye, Lay, et al. (2016) for planar megathrust faults with uniform dip angles ranging from 21° to 23° using teleseismic body waves, we found slip patches near 20 and 10 km depth, respectively. Given that uniform dip extending to very shallow depth is unlikely, we conduct inversions of teleseismic body waves for a 2.5D fault surface, which varies in dip with depth along the fault, guided by the SLAB2 plate interface model (Hayes et al., 2018), along with variable bathymetry above the fault (Ye, Kanamori, et al., 2016), as there are only minor variations along strike over ∼100 km. This model has a strike of 276.9°, and the dip increases from 5° near the trench to 37° near 30 km depth (Figs. 2 and 3). The water depth varies from 6.0 km near the trench to 1.1 km above the down‐dip edge of the model (Fig. 3).

The seismic waveform data used for the inversion are 70 P waves and 33 SH waves from global broadband seismic stations in the distance range of 30°–95° with good azimuthal distribution (Figs. 2d and 4). The raw seismic waveform recordings are corrected by removing instrument response to ground displacement bandpass filtered in the frequency band of 0.005–0.9 Hz. We align the traces on very weak initial arrivals consistent with the USGS‐NEIC origin, with there being about 10 s before strong arrivals that are captured in the inversion.

We apply a least‐squares kinematic inversion with multiple rake‐varying (±45° from a reference value given by our W‐phase inversion) subfault source time function windows (e.g., Hartzell and Heaton, 1983; Kikuchi and Kanamori, 1991; Ye, Lay, et al., 2016) to determine the space‐time slip distribution for the $Mw$ 7.7 event using the 2.5D dip‐varying fault grid (Figs. 2 and 3). Model Crust 1.0 (Laske et al., 2013) was used for the local source velocity structure, adapted to have variable bathymetry over the fault. Local 1D Green’s functions are used for each row of the fault, approximating the far‐field P and SH response for the 2.5D structure.

For our preferred model, the grid has 22 7.5 km × 7.5 km subfault columns along strike and 10 rows along dip, with the subfault source time functions being represented by 12 2.0 s rise time‐symmetric triangles, with 2.0 s time shift between triangles, giving maximum possible subfault duration of 26 s. An upper bound of 2.5 km/s was imposed for the rupture speed. A hypocenter at 22.93° S, 171.7° E, 8.8 km depth was assumed to reconcile the USGS‐NEIC epicenter with the 2.5D fault grid (blue star in Fig. 3). The result of the 2.5D inversion of seismic data is shown in Figure 2, with seismic moment, $M0=4.17×1020N·m$ (⁠$Mw$ 7.68), and a slip centroid depth of 17.6 km. Peak slip of about 4.2 m is found in a concentrated slip patch that extends bilaterally over ∼100 km along strike. The moment rate function commences about 12 s after the origin time, with unresolved early slip, and then has about 50 s duration, with two peaks at around 25 and 35 s into the rupture originating from the bilateral growth of the slip zone. We find almost no slip at depths shallower than 10 km, when we used the more realistic interplate geometry. The radiated energy $ER$ estimated from the moment‐rate function for periods longer than 20 s and from the average of propagation‐corrected broadband P waves for periods from 20 to 1 s is $2.38×1015J$⁠, giving a moment‐scaled radiated energy (⁠$ER/M0$⁠) of $5.7×10−6$⁠, somewhat low among interplate thrust events (average $ER/M0∼1.1×10−5$⁠), but not as low as observed for tsunami earthquakes with $ER/M0$ ranging from $∼1×10−6$ to $4×10−6$ (Ye, Lay, et al., 2016).

The waveform fits are very good, with the model accounting for about 86% of the power in the observed signals (Fig. 4). The waveforms are relatively smooth and the two‐subevent character of the moment rate function is reflected in the double arrivals evident in the P waves at various azimuths. Comparable waveform fits and similar slip distribution are found for models with maximum rupture speed of 2.0 and 3.0 km/s, so the resolution of the bilateral rupture expansion speed is limited, but the basic depth range and slip extent are stably determined. The static stress‐drop estimates for this model are 5.5 MPa for a slip‐weighted estimate and 5.3 MPa for an equivalent circular model with area corresponding to the subfaults with more than 15% of the peak slip following Ye, Lay, et al. (2016).

Large tsunamis have seldom been generated by earthquakes in southern Vanuatu. The 28 March 1875 M ∼8.1 and 1920 M ∼7.8 events located well to the northwest from the 2021 sequence are among the most damaging (Sahal et al., 2010; Ioualalen et al., 2017; Fig. 1a). The 2021 $Mw$ 7.7 mainshock produced modest tsunami recorded at distant NOAA DART stations, tide gauges stations in New Caledonia and New Zealand, and New Zealand DART stations along Tonga and the southern New Hebrides. The tsunami waves at these stations are mostly less than 20 cm in amplitude, consistent with expectations given the moderate size of an event with slip of less than 4.2 m at depths from 10 to 25 km.

The up‐dip extent of slip is always difficult to resolve, and is important for the strength of tsunami excitation. If significant slip extends to shallow depths close to the trench, strong water multiples (pwP) will be excited in the deep water, giving ringing P‐wave coda that will not be as strongly excited if slip is confined to greater depth (e.g., Lay et al., 2019). For the $Mw$ 7.7 mainshock, we examined azimuthally distributed P‐wave coda amplitudes relative to direct P amplitudes at epicentral distances from 80° to 120° and found elevated median ratios (∼0.79) relative to events for which no slip occurred on the shallow megathrust but with slip typically in the 20–40 km depth range (fig. 6 in Lay et al., 2019). The median $mB$ of 7.27 for a period of 4.9 s (51 channels from 30° to 80°), and the coda magnitude $mcoda$ median value of 7.11, give $mcoda−mB=−0.16$⁠, which is within the range of values for events that do have some shallow slip in figure 11 of Lay et al. (2019). These observations are consistent with the primary slip locating relatively shallow on the megathrust, but not as concentrated below deep water as for the uniform fault dip inversions. In the depth‐varying megathrust scenario of Lay et al. (2012), this earthquake appears to primarily rupture Domain B (15–35 km deep), with minor penetration up into, but not total rupture, Domain A (<15 km deep).

The USGS‐NEIC distribution of seismicity prior to and after the large 2021 southern Loyalty Islands event is shown in Figure 1a. The aftershock sequence, involving 347 aftershocks with $Mw≥4.5$ reported by USGS‐NEIC in the first 45 days, is highly productive. Globally, shallow earthquakes with $Mw$ 7.7 on average have about 62 events in that magnitude and time interval (Ye et al., 2020), and this sequence is at the upper limit of the global observed aftershock productivity distribution. Four of the aftershocks have $Mw$ in the range 6.1–6.3, with two being thrusting events and two, including the largest event, having normal‐faulting mechanisms. One hundred and nine aftershocks in the first 45 days have $Mw≥5.0$⁠. Many aftershocks are located south of the trench with the larger events involving shallow normal faulting (Figs. 1 and 4). The region of large‐slip zone for the 2021 rupture has only a few overlapping Global CMT solutions from 1976 to 2020 and is located down‐dip of a band of shallow thrusting mechanisms, although uncertainty in the centroid locations must be kept in mind. Prior seismicity from the USGS‐NEIC catalog has shallow activity as well, but there is some overlap with the up‐dip portion of the 2021 large‐slip zone.

Strong activation of outer trench slope tensional ruptures in the aftershock sequence of an interplate thrust event is closely associated with shallow slip occurring on the megathrust (e.g., Wetzler et al., 2017; Sladen and Trevisan, 2018; Ye et al., 2021). Interaction of the stress regime in the outer rise and trench slope with the interplate earthquake cycle has long been discussed (e.g., Christensen and Ruff, 1988), and recent updates strongly support the increase in shallow tensional faulting after large interplate thrust earthquakes overall (Ye et al., 2021). The basic explanation is that when the interplate boundary is locked and accumulating strain, the plate bending stress near the trench is reduced, suppressing tensional faulting. Release of the strain accumulation with megathrust rupture abruptly allows elastic bending stresses to increase, producing normal faulting seaward of the large‐slip zone in the megathrust rupture. This model is generally consistent with the observations for the 2021 $Mw$ 7.7 southeast Loyalty Islands sequence.

Strong interactions between the megathrust and intraplate environment seaward of the trench have also been observed for numerous sequences along the southern Vanuatu arc and the Loyalty Islands, northwest of the 2021 sequence. Figure 5 shows time series of the 1995, 2003–2004, and 2018 sequences, including both USGS‐NEIC seismicity and Global CMT focal mechanism catalogs. They all have significant aftershock counts within the magnitude‐scaled source region (white circles) for $Mw≥4.5$ events within 45 days of the mainshock, suggesting that these events are all more productive than the global averages for similar size mainshocks (Ye et al., 2020). Together with the 2021 sequence, the southern Loyalty Islands region is found to have characteristically high aftershock productivity and unusually strong interactions between megathrust and outer trench slope activity in general.

The 2021 southeast Loyalty Islands interplate thrust sequence involved a very productive faulting activity with large foreshocks and many aftershocks of the $Mw$ 7.7 mainshock. Slip during the mainshock was concentrated at relatively shallow depths of 10–25 km, but it does not appear that much slip occurred in the uppermost portion of the megathrust. The shallow slip generated a modest tsunami detected as far away as New Zealand, but the rupture does not appear to be as slow or as energy deficient as would characterize a tsunami earthquake. The relatively shallow concentrated slip does appear to have enhanced triggering of outer trench slope normal faulting directly seaward of the slip zone, producing several aftershocks larger than $Mw$ 6.0 and many smaller events. This interaction between the megathrust and the intraplate environment seaward of the trench is a robust characteristic of the productive earthquake sequences along the Loyalty Islands.

Body wave and surface wave recordings from global seismic stations were downloaded from the Incorporated Research Institutions for Seismology (IRIS) Data Management Center (https://ds.iris.edu/wilber3/find_event, last accessed November 2021). Global Centroid Moment Tensor (Global CMT) solutions are from https://www.globalcmt.org/CMTsearch.html (last accessed November 2021). Earthquake information is based on the catalog of the National Earthquake Information Center at the U.S. Geological Survey (USGS‐NEIC) (https://earthquake.usgs.gov/earthquakes, last accessed May 2021). The hypocenter and finite‐fault information for the 10 February 2021 $Mw$ 7.7 mainshock from the USGS‐NEIC are available at https://earthquake.usgs.gov/earthquakes/eventpage/us6000dg77/executive (last accessed November 2021).

All data used in this study are openly available from the sources listed in Data and Resources.

The authors acknowledge there are no conflicts of interest recorded.

The authors thank Editor Keith Koper, Associate Editor Ruth Harris, and two anonymous reviewers for their helpful comments on this article. Lingling Ye’s and Xiaofei Chen’s earthquake study are supported by National Natural Science Foundation of China (41874056, 41790465, and U1901602). Thorne Lay’s research on earthquake processes is supported by National Science Foundation Grant EAR1802364.