Relocation of the April 2017 Batangas, Philippines, earthquake sequence, with tectonic implications

The 2017 Mw 5.9 Batangas earthquake occurred on a location abutting the Macolod Corridor (MC) to the NE and above the steeply dipping slab of subducted South China Sea lithosphere to the SW. While the spatial distribution of the 2017 Batangas earthquake sequence provides a constraint on the scope of MC’s geothermal activity, the focal mechanisms are associated with the states of stress in MC. Here, we relocated the 2017 Batangas sequence using data from two seismic network which cover the region of the NW Mindoro and Southern Luzon. Upon employing the HypoDD algorithm, we retrieve a compact distribution within the Batangas Bay, which in turn suggests that the Balayan Bay to the north remains aseismic and is within the scope of geothermally impact area by MC. In addition, the 2017 Batangas earthquake and the 1994 Mindoro Earthquake share similar focal mechanisms that exhibit an E-W extension in the region of SW Luzon, which differs from a few leftlateral moderate earthquakes to the east at 121.5°E on the Sibuyan Verde Passage Fault. Together with the results of 3D numerical modeling on coeval subduction of oceanic and continental lithospheres and P-wave tomography in the region from previous studies, we incline to explain the regional E-W extension in SW Luzon as induced by the toroidal asthenospheric flow and the mechanisms of MC formation as originated by the overriding-plate stretching. Both are incipiently driven by the collision of the Palawan Continental Block into the Philippine Mobile Belt and the subsequent detachment of the South China Sea slab beneath central Mindoro. Article history: Received 20 September 2019 Revised 30 January 2020 Accepted 31 January 2020


INTRODUCTION
The 8 April 2017, M w 5.9 Batangas earthquake at 16.1 km centroid depth was preceded by a prominent foreshock sequence starting with a M w 5.2 earthquake on 4 April as recorded by the Global Centroid Moment Tensor (GCMT) catalog (Dziewonski et al. 1981). The largest earthquake in the region during the past few decades was the 14 November 1994, M w 7.1 Mindoro tsunamigenic earthquake (Tanioka and Satake 1996) (Fig. 1a), which hypocenter was located at (121.1°E, 13.5°N) by the Philippine Institute of Volcanology and Seismology (PHIVOLCS) and at (121.067°E, 13.525°N) by NEIC Preliminary Determination of Epicenters (PDE) and which centroid was located about 20 km to the east at (121.32°E, 13.44°N) in GCMT catalogue. The surface rupture of the 1994 Mindoro earthquake has been mapped along a previously unmapped, N-S Aglubang River Fault (ARF) with dominantly right-lateral strike slip over a length of about 35 km in NE Mindoro (PHIVOLCS Quick Response Teams 1994). Main shocks of the 2017 Batangas and 1994 Mindoro earthquake share similar focal mechanisms with a nearly E-W tension axis (Fig. 1a). The major fault in the region is the WNW-ESE Sibuyan Verde Passage Fault (SVPF), which branches out from the Philippine Faults (PF) near Masbate Island and passes between Mindoro Island and Luzon westward to the Manila Trench (Marchadier and Rangin 1990) (Fig. 1b). Like PF, the SVPF is broadly in left-lateral motion, absorbing the oblique collision between the Philippine Sea plate and the Sundaland block (Galgana et al. 2007) as attested by a few ~M w 5.8 left-lateral earthquakes to the east at 121.5°E around SVPF. However, if the WNW-ESE trending SVPF was taken as the fault plane for the 2017 Batangas and 1994 Mindoro earthquakes, these fault motions would be right-lateral movements, suggesting the necessity of local E-W tensional stress different from those of the SVPF in a regional scale.
Location of the 2017 Batangas earthquake also bears significant tectonic implications. To the northeast, the 40 km-wide NE-SW Macolod Corridor (MC) is a transtensional relay zone expressed as left-lateral pull-part grabens and active volcanism (Förster et al. 1990). The fact that the MC is currently under E-W tensional stress but devoid of shallow seismicity (Fig. 1a) is attributed to a hot geotherm caused by the Taal/Laguna de Bay volcanic field (Pubellier et al. 2000). The high temperature may weaken the rigidity of the lithosphere, making the area in the vicinity of MC incapable of accumulating tectonic strain for earthquake ruptures. In this context, the occurrence and distribution of the 2017 Batangas earthquake sequence highlight the area where the crustal deformation is beyond the geothermal impact of MC and, thus in turn, constrain the spatial scope of MC's geothermal activity. In addition, the proximity of the 2017 Batangas earthquake and MC would be exposed under the same regional stress patterns, which can be learnt from the focal mechanisms of the main shock. To the southwest, on the other hand, the 2017 Batangas earthquake sits almost on top of a steeply dipping slab from the subduction of the South China Sea lithosphere along the Manila Trench (Chen et al. 2015;Bina et al. 2020) (Figs. 1a and c). Here, the Manila Trench was terminated southward at Mindoro Island by the indentation of the Palawan Continental Block (PCB) into the Philippine Mobile Belt (PMB) during the Miocene (Yumul et al. 2003;Doo et al. 2015). The geodynamic evolution of coeval subduction of oceanic and continental lithospheres was recently modeled by 3D numerical simulations to investigate magma genesis and transport processes associated with slab roll-back and tearing (Menant et al. 2016). Together with results of 3D P-wave tomography in the region, we provide an explanation on the regional E-W tensional stress in SW Luzon and mechanisms of MC formation in discussion.
Since 2013, we have deployed five broadband stations in NW Mindoro (marked by red triangles in Fig. 2) to monitor seismic activity in the region and data from the 2017 Batangas earthquake sequence were well recorded. The additional data significantly reduce the azimuthal gap of stations and alleviate the tradeoff between distance and depth caused by events outside of the PHIVOLCS network. In combination with the data from six nearby PHIVOLC's stations (blue triangles in Fig. 2), we employed the doubledifference location algorithm (HypoDD; Waldhauser and Ellsworth 2000) to simultaneously relocating events in the same cluster and together to determine the distribution of the Batangas sequence. The seismicity rate and magnitudefrequency relationship of both foreshock and main shock sequence were then analyzed. We aim at understanding the tectonic implications and seismicity characteristics of the 2017 Batangas earthquake sequence through a thorough study of its spatial and temporal distributions.

DATA AND METHODS
To auto-pick seismic signals in waveform data from the combined seismic network (triangles in Fig. 2), we invoked the algorithms of short-term average over long-term average (Allen 1982). Having bandpass filtered the data between 1 April and 22 April 2017 with a frequency band between 1 and 10 Hz, the short-term and long-term moving averages (STA and LTA) of waveform amplitude were calculated using a 0.5-s and a 10-s time window, respectively. A seismic signal was picked once STA/LTA ratio exceeds 4 and keep at least 2 for more than 2 s. These auto-picked arrivals were assigned to trial hypocenters by grid-searching a region parameterized by grid points with intervals of 0.01° in both N-S and E-W direction and 2 km from 0 to 50 km depth and 5 km from 50 to 100 km depth in depth direction and within the following boundaries: 119.92 -122.92°E, 12.68 -14.98°N, and 0 -100 km depth. For each event, stations with either P residuals less than 1.5 s or S residuals less than 3.0 s were kept. Only those events with at least four kept stations within 1° or six stations within 2° epicentral distances were kept for further analysis, leaving a total of ~3000 events. For these events, we visually inspected their P-and S-waves in the time window of predicted arrivals of all stations and retained those events providing clear onsets of P-and S-waves, leaving a total of 361 events. We manually re-picked their arrivals with assigned weightings and located them by implementing an iterative least-square method (Pavlis et al. 2004) through a graphical user interface (dbloc2) in Antelope Software Package (http://www. brtt.com). The 1-D velocity structure ( Fig. 2) used to locate the earthquakes was determined by using data recorded by short-period network deployed on NW Mindoro and SW Luzon during April 2010 to April 2011 (squares in Fig. 2). 169 events (Fig. 2) located by using at least four P readings and one S reading from the network were selected to invert the 1D velocity structure by applying the VELEST method (Kissling et al. 1994).
Next, the 361 events of the 2017 Batangas earthquake sequence were relocated by applying the HypoDD algorithm, which iteratively invert for hypocenter locations of closely spaced event pairs by minimizing the residuals of their travel-time differences. The procedures effectively minimized the location errors induced by velocity structures. We first formed event pairs whenever neighboring events within 5 km had at least 8 common observations; we kept only the 15 nearest event pairs for each event. Then, the residual differences of common observations for each event pair were calculated. A total of 3342 event pairs with 19881  P-wave and 18240 S-wave residual differences were obtained from the 361 events. Lastly, a weighted least-square inversion was conducted for 20 iterations to determine the optimal hypocenter locations. The damping factor dropped from 50 to 45 after the 10 th iteration and the weighting was set to zero for outlying observations after the 5 th iteration. The root mean square of the travel-time difference residuals reduces from 0.379 to 0.207 s after HypoDD relocation.
Having relocated earthquake hypocenters, we determined the local magnitude M L using the maximum amplitude A (in μm) and epicentral distance D (in km), assuming a linear relationship M L = log 10 A + c 1 + c 2 [log 10 D] (Lay and Wallace 1995). The intercept (c 1 ) and slope (c 2 ) were calibrated using three stations' data from three events of known M w that generally exhibit a linear trend (Fig. 3). This scaling of M L with M w allowed for M L -moment conversion using log 10 M 0 = 1.5M w + 9.1, based on the Gutenberg-Richter magnitude-energy relation (Gutenberg and Richter 1944) and an approximate relation between radiated wave energy and seismic moment (Kanamori 1977).

RESULTS
Comparing distributions of seismicity determined by the HypoDD algorithm with those located by PHIVOLCS, the former clearly shows more clustering than the latter (Fig. 4). Also, hypocenter location as downloaded from the PHIVOLCS website has lower precision, 0.01° and 1 km in horizontal and vertical direction, respectively, as shown by the grid-like distribution of the seismicity (Fig. 4). Precision of origin time from PHIVOLCS is 1 s. Using events in the PHIVOLCS catalog as reference and searching the HypoDD catalog within ±30 s of origin time and 50 km from the hypocenters, the geographically nearest HypoDD entry was taken as the same event and connected in order to investigate shifting of hypocenters between the two catalogs. Results exhibit a clear trend that most events in the PHIVOLCS catalog moved southward and upward to their counterparts in the HypoDD catalog (Figs. 4 and 5). The trend can be attributed to the fact that addition of Mindoro stations to the south makes events within the network and helps constrain depth. We illustrate this concept by drawing a circle with PHIVOLCS station (TGYB) at center and a proper radius to follow the trend (Fig. 4b). For events locating on the circle, the hypocentral distances are all the same for stations near the center. In other words, one-sided stations suffer from tradeoff between distance and depth, which can be alleviated by adding stations on the other side where the hypocentral distances are different on the circle. To further show that the trend is mainly due to the addition of Mindoro stations, we performed an iterative leastsquare method on data of both the PHIVOLCS and Mindoro networks. Result already show a concentration within the Batangas Bay. Further application of HypoDD only made a minor improvement (Supplementary, Fig. S1), suggesting the vital role of Mindoro stations on earthquake relocations. Source parameters of the HypoDD catalog are listed in a supplementary table (Supplementary, Table S1).
Next, we examined the temporal and spatial distributions of the Batangas sequence in the HypoDD catalog. Events of M L ≥ 5.0 or M w ≥ 5.0 were first ordered and marked chronologically. Second, the events in the foreshock sequence were displayed separately from those of the main shock sequence in Figs. 6a and b, not only in map view but also in cross sections projected on two nodal planes of the largest foreshock (No. 1) and the main shock (No. 7), respectively. Both sequences elongate in a direction favoring the NNW-SSE nodal plane as the fault plane. The foreshock sequence is distributed along ~8 km of length and ~18 ± 5 km of depth, whereas the main shock sequence is distributed along ~20 km of length and ~17 ± 6 km of depth. While the similar depth ranges may be controlled by the seismogenic depth, the length difference is a reflection of earthquake size difference (M w 5.2 versus M w 5.9). It appears that most M L ≥ 5.0 events occurred at deeper depth, except for the main shock. As for temporal evolution of seismicity, the cumulative moment release as a function of time was calculated using a bin size of 3 hours (Fig. 7). Results show that the foreshock sequence initiated with intensive events in the first 9 hours and then gradually decayed for about three days until the occurrence of the main shock.

DISCUSSION
Incorporation of the data from broadband stations in NW Mindoro fills the station network gap and alleviates the tradeoff between distance and depth caused by events outside of the network. As a result, most events originally located in Balayan Bay at greater depths in the PHIVOLCS catalog were moved to Batangas Bay at shallower depth in the HypoDD catalog, forming a more compact cluster. The picture better approximates the typical ~100 km 2 empirical fault area for a M w 6.0 event (Wells and Coppersmith 1994). More importantly, the refined relocations of earthquake events were not distributed beyond Batangas Bay. This in turn suggests that Balayan Bay to the north remains largely free of shallow seismic events and is within the scope of the geothermally impact area. A high resolution velocity structures may help understand the seismogenic structures in the region (Wen et al. 2019), which needs more dense station coverages.
Whether the NNW-SSE trending fault of the 2017 Batangas earthquake represents the northern extension of ARF of the 1994 Mindoro earthquake is unknown. Focal mechanisms of both earthquakes exhibit E-W extension, as well as the E-W rifting of MC, suggesting the dominant E-W tensional stress patterns in SW Luzon. Here, previous dynamic studies on simulations of continental subduction/ Only data that follow the linear trend were used (three events at three stations). The event numbers (6,7,8) can be used to find the event locations in Fig. 6 and the orders (102, 103, 106) for table in Table S1.  collision were invoked to explain stress patterns of SW Luzon and mechanism of MC. The 2D modeling on scenarios of the Mindoro Island (Bina et al. 2020) reveals that, upon incipient collision of approaching continents, trench retreat is succeeded by trench advance and slab dip begins to increase due to the combined effects of the negative thermal and petrological buoyancy, which eventually induce slab weakening, necking, and detachment. Putting into 3D scenarios for coeval subduction of oceanic and continental lithospheres (Menant et al. 2016), the along-trench buoyancy variations and the slab tear under the continental collision would induce stretching of the overriding plate above the oceanic subduction and a toroidal asthenospheric flow. Recent P-wave tomography indeed exhibits a detached slab and apparent tearing in central and north Mindoro, respectively (Fan et al. 2017). Using high Vp anomalies for subducting slab ) and low Vp anomalies for positions of slab detachment beneath Mindoro and SW Luzon (Fan et al. 2017) and invoking results of 3D numerical modeling (Menant et al. 2016), we were able to infer the pattern of toroidal asthenospheric flow and location of overriding-plate stretching in the study region (Fig. 8). While the flow pattern may induce E-W tensional stress in SW Luzon, the stretching may be responsible for MC formation. We thus conclude that the regional E-W extension in SW Luzon and the mechanisms of MC may incipiently be driven by the PCB-PMB collision.
In terms of seismicity characteristics, the 2017 Batangas earthquake exhibits prominent foreshocks immediately before the main shock. The feasibility of using immediate foreshocks to forecast the ensuing main shock remains an open question, as their causal and statistical relationships remain unclear (Dodge et al. 1995;Felzer et al. 2004). One of the often investigated potential precursor is the difference in b-values between foreshock and aftershock sequences (Knopoff et al. 1982). We conducted the magnitude-frequency analysis for both the foreshock and aftershock sequences to determine the b-value by linear fitting (Fig. 9), where the magnitude of completeness was derived by the point with the maximum value of the first derivative of the frequency-magnitude curve (maximum curvature; Wiemer and Wyss 2000). The b-values thus derived were   comparable with those determined by the maximum likelihood method (Woessner and Wiemer 2005;Chen et al. 2008). We found that the b-value is ~0.6 for the foreshock and ~0.7 for the aftershock sequences (Fig. 9). The drop of b-value for the foreshock sequence was also observed for the 6 February 2018, M w 6.4 Hualien earthquake Rau and Tseng 2019), suggesting a high seismicity rate with events of relatively large magnitude and may be potential for forecasting the issuing main shock (Gulia and Wiemer 2019). The b-values less than typical value (one) could due to the nature of foreshock and aftershock sequences as manifested by most aftershock sequences in east Taiwan  or due to the general low b-values in SW Luzon as revealed by the comprehensive survey of the entire Philippines archipelago (Pailoplee and Boonchaluay 2016).

CONCLUSIONS
We relocated the earthquakes of the 2017 Batangas sequence using additional data from five broadband stations in NW Mindoro and employing the HypoDD algorithm. The resulting distributions of seismicity exhibit a more compact NNW-SSE trending cluster limited within Batangas Bay. The focal mechanism of the main shock appears similar to that of the 1994 M w 7.1 Mindoro earthquake, manifesting E-W extension in the region of SW Luzon, which also controls the rifting direction of the MC. Based on the results of 3D numerical modeling on coeval subduction of oceanic and continental lithospheres and of 3D P-wave tomography in the region from previous studies, we propose that regional tensional stresses are likely induced by the toroidal asthenospheric flow and mechanisms of MC formation are likely resulted from overriding-plate stretching. Both evolved from slab tearing due to subduction/collision of the PCB into PMB. Lastly, seismicity characteristics of the 2017 Batangas earthquake exhibit a prominent foreshock sequence initiated about four days before the main shock, with a b-value smaller than that of the aftershock sequence.