Source Model and Characteristics of the 27 July 2022 M W 7.0 Northwestern Luzon Earthquake, Philippines

Thegeometryandkinematicsofthecausativefaultofthe27July2022momentmagni-tude (M W ) 7.0 earthquake, which is one of the strongest to hit northern and central Luzon in the past 30 years, were estimated through inverse modeling of line-of-sight interferometric synthetic aperture radar deformation. We modeled rupture along multiple candidate faults based on fit with the pattern of line-of-sight deformation, consistency with focal mechanisms, and compatibility with the known kinematics of the mapped active faults in the region. Our preferred fault model, located west of and parallel to the Abra River Fault (ARF), exhibits localized reverse-slip (average 67° rake) at 15-35 km down-dip. Peak slip occurs at 13-16 km depth, with 95 cm of pure reverse-slip. The existence of a reverse-slip dominated ARF-parallel fault rupture is consistent with a complex shear partitioning model, wherein the NW-SE oblique plate convergence is accommodated not only by the sinistral strike-slip Philippine Fault Zone and the major subduction zones, but also by minor faults in intervening crustal blocks.


Introduction
On 27 July 2022 at 08:43 local time (UTC +8), the northwestern region of the island of Luzon, northern Philippines was hit by a moment magnitude (M W ) 7.0 earthquake ( Fig. 1 A and B). The epicenter was located 10 km south of Tayum, Abra (17.5°N, 120.7°E), and had a focal depth of~20 km ( Fig. 1 B); (PHIVOLCS, 2022a,b). Focal mechanisms produced both by the Philippine Institute of Volcanology and Seismology (PHIVOLCS) and the United States Geological Survey (USGS), consistently suggest oblique-reverse faulting on either a N-striking, E-dipping or a SW-striking, NW-dipping fault ( Fig. 1 B); (PHIVOLCS, 2022b;USGS, 2022). The PHIVOLCS focal mechanism has a strike of 8°, a dip of 28°, and rake of for this event has yet to be located and properly mapped.
A total of 11 fatalities, 609 injuries, 49,803 displaced persons, and~US$45 million in damages to major public infrastructure and to the agricultural sector was reported (NDRRMC, 2022). The intensity values and distributions, though slightly different depending on the reporting agency, are in agreement with the highest felt intensity being centered in the province of Abra. A modified Mercalli intensity value of 7.5 was reported by the USGS, and a Philippine earthquake intensity scale (PEIS) value of 7 was reported by PHIVOLCS, which also included the western coastal towns in the province of Ilocos Sur (USGS, 2022;PHIVOLCS, 2022d).
In the past 50 years, 11 M W > 6.5 earthquakes have occurred within 250 km of the 2022 Luzon epicenter in Tayum, Abra (USGS, 2022). The largest earthquake in history to occur in northern Luzon was the 1990 M W 7.7 Luzon Earthquake, which was associated with a~120 km-long surface rupture along the Digdig Segment of the Philippine Fault Zone (Punongbayan et al., 1991-07-16;Nakata et al., 1996).
This study presents the first interferometric synthetic-aperture radar (InSAR)-based fault source model for the 2022 M W 7.0 northwestern Luzon earthquake, providing information such as the possible location, geometry, and slip distribution. Such data can contribute to a better understanding of this particular event, as well as generally the styles, mechanisms, and distribution of deformation in the Philippines-a tectonically complex, seismically active region which could benefit from a more comprehensive mapping and accurate kinematic analysis of active structures. This kind of effort would enhance the country's capability to assess seismic hazards and risks.

Tectonic Setting
The 8.0 cm/yr northwestward motion of the Philippine Sea Plate (PSP) towards the Sunda Plate (SP) (Seno et al., 1993) is accommodated throughout the Philippine archipelago by a system of crustal faults and subduction zones that exhibit complex shear partitioning ( Fig. 1 A; e.g., Rimando et al. (2019). To its west and east, the island arc is bound by the east-dipping Manila-Negros-Sulu-Cotabato Trench System and the west-dipping East Luzon Trough-Philippine Trench System, respectively. In between these trenches, is the 1400 km-long, sinistral Philippine Fault Zone (PFZ), which runs along the entire length of the archipelago, from the island of Luzon in the northwest to the island of Mindanao in the southeast (Allen, 1962;Hamilton, 1979;Acharya and Aggarwal, 1980;Bautista et al., 2001;Cardwell et al., 1980;Hamburger et al., 1983;Hayes and Lewis, 1985;Ozawa et al., 2004;Rimando and Knuepfer, 2006;Marfito et al., 2022) (Fig. 1 A). There is an estimated 80-100 km and 200 km of minimum displacement along the PFZ in northwest Luzon (Pinet and Stephan, 1990) and Mindanao (Mitchell et al., 1986), respectively, since the Miocene.
The boundary-perpendicular component of the overall oblique plate convergence is accommodated by subduction zones, inferred thrust/reverse and oblique strike-slip faults in the crustal blocks bounded by major active faults, and by regional tectonic uplift, while the boundary-parallel component is accommodated mostly by the PFZ (Fig.1 A; e.g., ). In northwestern Luzon, however, the Vigan-Aggao Fault (VAF), which forms the westernmost strand of the PFZ, also accommodates a significant portion of trench-perpendicular shortening through angled sinistral strike-slip faulting (e.g.,  and Fig. 1 B). While other minor active faults likely exist within the crustal blocks, the exact traces and kinematics of these have yet to be comprehensively documented.

InSAR processing
We used the descending track (20220721-20220802) Sentinel-1A synthetic-aperture radar (SAR) single-look complex images from the European Union's Copernicus Programme satellite constellation to create an interferogram and a line-of-sight (LOS) displacement map of the area between six days before and six days after the event. Unfortunately, there were no acquisitions for the ascending track, due to the end of Sentinel-1B's mission in late 2021 (ESA, 2022). Sentinel-1A uses the C-band, corresponding to a wavelength of 5.5 cm. The images were acquired in 'Terrain Observation with Progressive Scans' mode, which bundles three sub-swaths together to cover a greater area. However, given the limited size of the earthquake and distribution of deformation, only the westernmost sub-swath was processed and analyzed.
We used the Generic Mapping Tools Synthetic Aperture Radar (GMTSAR), an open source InSAR processing program (Sandwell et al., 2011), to carry out our analysis. The SAR images and their precise orbital information were obtained from the Copernicus Open Access Hub (Copernicus, 2022). A 30 m-resolution Shuttle Radar Topographic Mission Version 3 (SRTM1v3) digital elevation model (DEM) of the area was generated using the online GMTSAR DEM Generator (G.M.T.S.A.R., 2010), and was used to correct for topography. A Gaussian filter with a wavelength of 200 m was applied to the images, and the pixels were decimated by a factor of two along the azimuth and by a factor of eight along the range prior to creating the wrapped interferogram (Sandwell et al., 2011). A coherence mask with a threshold of 0.085 was applied to the data prior to unwrapping with the SNAPHU algorithm (Chen andZebker, 2000, 2001;Chen, 2002). The unwrapped phase was corrected for tropospheric effects using data from the General Atmospheric Correction Online Service for InSAR (Yu et al., 2017(Yu et al., , 2018a before converting to LOS displacement (Fig. 2).

Earthquake Source Modeling
We solved for earthquake rupture using two simple planar fault geometries ( Fig. 3 A and B), the selection of which is guided by the PHIVOLCS and USGS focal mechanisms, the visible pattern of line-of-sight deformation from the unwrapped interferogram, and the mapped active faults in the region. Based on these considera-tions, we explored the possibility of rupture along and parallel to two local faults: the Vigan-Aggao Fault (VAF) and the Abra River Fault (ARF), respectively. There are other faults that are local to the rupture, includ- ing the Naglibacan Fault and Bangui Fault; however due to their sub-optimal orientations in relation to the PHIVOLCS focal mechanism, these were not considered to be ideal candidates for rupture and therefore were not tested. The well-documented NNE-SSW Vigan-Aggao Fault (VAF), is a range-bounding fault that is located close to the western coast of northern Luzon, and parallels both the trend of the hinge between the positive LOS and negative LOS deformation, and the trend of the aftershock distribution. This active fault is known to have primarily left-lateral displacement, and has been active since the Pliocene (Pinet and Stephan, 1990). Although the Luzon earthquake of 2022 was primarily a thrust mechanism, we allow for the possibility of rupture on the VAF due to its proximity. We also selected a more N-S fault plane to parallel the more proximal but less well studied Abra River Fault (ARF) (Pinet and Stephan, 1990). Since the ARF has been described as a dominantly strike-slip fault (Pinet and Stephan, 1990), rupture along the ARF itself, which has been suggested by PHIVOLCS to be the likely causative fault (PHIVOLCS, 2022d), is unlikely. Therefore, we modeled a fault plane, which fits the LOS deformation, to the west of the ARF. While an investigation of the optimum fault plane was beyond the scope of this rapid study, future efforts may be able to better constrain the ideal fault plane for this event based on geodetic data and any additional information about surface rupture, if it is found. Both faults are modeled with 30°E dips, consistent with the USGS and PHIVOLCS focal mechanism solutions. For comparison, we also used the fault plane solution provided by the USGS (USGS (2022); Fig. 3 C), and modeled the expected LOS displacement at the same locations as used in our source inversion (Fig. 3 D, E and F). The USGS finite fault model uses teleseismic body and surface waves and follows the methods of Ji et al. (2002). The parameters of all three faults are presented in Table 1.
To invert for the slip, we used the MudPy modeling and source inversion toolkit (Melgar and Bock, 2015). Because InSAR data is insensitive to the earthquake rupture velocity, we solved for slip as a static rupture. Green's functions for the InSAR data are calculated using the frequency wavenumber methods from Zhu and Rivera (2002). We used a velocity model that is local to the epicenter through CRUST1 (Laske et al., 2013-04). Due to the tectonics of the region, left-lateral and thrust fault slip was enforced, limiting model slip to rakes in a window between 0°and 90°. The inversion results were constrained using a Tikhonov spatial regularization scheme (Mair, 1994;Tikhonov, 1963). This regularization scheme imposes equal amounts of smoothing across all subfaults in our model and is guided by a spatial regularization constant. As the constant approaches zero, the problem approaches a non-regularized least squares solution. We test our inversion results over a range of values, opting for the model that minimizes the data misfit without overfitting or allowing for too rough of a final solution (Fig.3 A and B). The preferred solutions sit close to the bend in an L-curve test (Fig. S1).

LOS Displacements
The LOS displacements, displayed in Fig. 2, are relative to a satellite heading of -167°looking right at an incidence angle of 44°. Results show a lobe of positive LOS deformation beneath the northern ARF, and a lobe of negative LOS deformation to the SSE. This satellite geometry is particularly well suited to image deformation in the dip-slip direction. Given focal mechanism results that suggest a N-striking fault plane, dipping to the E, the LOS displacements therefore suggest dominantly reverse motion, moving the eastern hanging wall vertically upwards.

Inversion Results
The inversion results for all three geometries tested are displayed in Table 1   dicted LOS displacements at the surface from a forward model of the slip at depth (Fig. 3, D-F), and the misfit between that forward model and the observed LOS displacements (Fig. 3, G-I

VAF rupture model
This model exhibits a diffuse amount of primarily reverse slip (average 63.6°rake) over an along-strike band between 32-40 km down-dip (16-20 km depth) with peak slip of 64 cm (Figs. 3A, S2 A, and S2 B). Lesser amounts of left-lateral slip are mainly concentrated on a narrow band along the central segment of the fault plane, starting at 60 km down-dip (30 km depth) and shallowing to the north (Fig. S2 C).

Abra River Fault-parallel model
This model exhibits more localized reverse-slip (average 67°rake) at shallower depths and with a higher peak slip value (95 cm) (Figs. 3 B, S3 A, and S3 B) than the VAF rupture geometry. Left-lateral slip is mostly confined to the central segment of the plane, at 20-40 km down-dip distance (10-20 km depth).

USGS forward model
The USGS forward model, which is based on the finite fault model that was released shortly after the main event, exhibits a very focused region of slip beneath the lobe of positive LOS displacement, with a more dominant left-lateral slip (average rake 34°) and with peak slip of around 80 cm occurring near 12 km depth (Figs. 3 C, F and I).

Discussion and Conclusions
Although both the VAF and ARF-parallel fault planes produced similar L2 misfits (Table 1) and are therefore similarly good choices, the modeled rupture on the ARF-parallel fault plane is the preferred model for this study as it visually most closely reproduces the overall observed LOS displacements near the expected earthquake rupture (Fig. 3 H). While the surface projection of this model is currently not associated with any mapped active fault trace, it is expected that there must be significant margin-normal shortening across the northern Philippines given the opposing subduction zones to the east and west (e.g., Bautista et al. (2001)). In this sense, such a fault is kinematically congruent with the VAF and ARF, which are both accommodating the mostly sinistral strike-slip component of the oblique convergence in northern Luzon. If this is the causative fault of the 2022 Mw 7.0 Luzon earthquake, there may well be surface rupture to the west of Abra, near the most productive aftershock region (Fig.1 B). It may also be worth considering the possibility of stress transfer onto the Abra River Fault, which would lie above the fault plane modeled herein.
On the other hand, rupture on the VAF, which is dominated by relatively deep reverse slip (peak slip between 30-40 km down-dip), is associated with the lowest misfit. It is worth noting, though, that there is a minimal difference in the misfit values of the VAF and ARFparallel fault models. Additionally, while the surface projection of the VAF model coincides with the trace of a well-known fault, the dominantly reverse-slip kinematics that our model suggests is contradictory to the dominantly strike-slip displacement that has been documented through detailed mapping and a quantitative analysis of morphotectonic kinematic indicators along this fault zone . If this fault were the causative fault for this earthquake, there may not be any surface rupture due to the paucity of strong slip in the shallow 0-25 km depth range.
Among the earthquake source models, the USGS forward model is associated with the highest overall misfit and strongest residuals (Fig. 3 I). The high misfit is due to the incompatibility of the modeled slip distribution with the observed positive LOS deformation that is visible further to the south. We therefore consider this as an unlikely causative fault for this event.
The likely presence of a fault parallel to and west of the ARF that exhibits dominantly reverse-slip further supports a complex shear partitioning model for the Philippine archipelago, wherein the oblique plate convergence in Luzon is accommodated not only by the sinistral strike-slip PFZ and the major trenches, but is also taken up substantially by faults within crustal blocks (e.g., ).
Additionally, the fact that a previously unmapped ARF-parallel fault may exist underscores the importance of a comprehensive mapping of active tectonic structures even in areas where there is very subtle or poor topographic expression of faulting. This could be achieved through a combination of high-resolution DEM inspection, InSAR analysis and modeling, highresolution subsurface imaging, high-resolution potential field surveys, seismotectonic analysis, measurement of geomorphic indices, and slip tendency analysis of geophysical or topographic lineaments.
In the absence of complete mapping of active seismotectonic structures, this earthquake serves as a reminder of the importance of comprehensive seismic hazard mapping that considers the effect of shallow crustal earthquakes on as-yet unknown faults. If this event were to have happened further to the south it may have caused much greater damages and losses, particularly if it were to have affected a major population center such as Baguio City.
Future studies may consider Bayesian inverse modeling that involves producing posterior probability distributions of source parameters of the fault rupture, for example through Markov-chain Monte Carlo-and Metropolis-Hastings-based algorithms of Bagnardi and Hooper (2018), to provide more constrained estimates of fault geometry and slip distributions. In the meantime, though, our models can serve as a guide for the ongoing search for surface rupture and/or nearby, preexisting yet unmapped potentially active faults. mica Fast Reports, and the anonymous reviewers. This work is funded in part by NSERC DG RGPIN-2021-03031, and is NRCan ESS contribution number 20220267. We wish the people of Abra province a quick and complete recovery from this event.
The interferogram data and slip model files are available in the open-access Zenodo repository: Source model inputs and results -The July 2022 Mw 7.0 Northwestern Luzon Earthquake, Philippines (http://doi.org/10.5281/zenodo.7033117).