The 2019-2020 Indios , Puerto Rico earthquake 5 sequence : seismicity and faulting

The 2019-2020 Indios, Puerto Rico earthquake sequence ruptured on multiple faults with 50 several moderate magnitude earthquakes. Here we investigate the seismotectonics of this fault system using high precision hypocenter relocation and inversion of the near-field strong motions of five largest events in the sequence (5.6≤Mw≤6.4) for kinematic rupture models. The Mw6.4 mainshock occurred on an NE striking, SE dipping normal fault. The rupture nucleated offshore ~15 km SE of Indios at the depth of 8.6 km and extended SW-NE and up55 dip with an average speed of 1.55 km/s, reaching the seafloor and shoreline after about 8 seconds. The Mw5.7 and Mw5.8 events on 6th and 7th January 2020, respectively, occurred on two E-SE striking, near-vertical, left-lateral strike-slip faults. However, the 7th January 2020 Mw5.8 aftershock which occurred only 10 minutes after the mainshock, ruptured on a fault with almost the same strike as the mainshock but situated ~8 km further E, forming a set of 60 parallel faults in the fault system. On 11th January 2020, a Mw6.0 earthquake occurred on a NNE striking, W dipping fault, orthogonal to the faults hosting the strike-slip earthquakes. We also apply template matching for the detection of missed, small magnitude earthquakes and to study the spatial evolution of the main part of the sequence. These detections show multiple examples of acceleration phases before moderate earthquakes and reveal migration patterns 65 within the sequence. Using the template matching results along with GPS analysis, we image the temporal evolution of a foreshock sequence (Caja swarm). We propose that the swarm and the main sequence were a response to an inferred tectonic transient that most likely originated on the Muertos subduction as a slow-slip event.


Introduction
The 2019-2020 Indios, Puerto Rico earthquake sequence included multiple moderate earthquakes (Mw>5.0) that ruptured different faults with different faulting mechanisms (University of Puerto Rico, 1986). The swarm-like Indios sequence started in July 2019 with a 75 smaller offshore Caja swarm, 14 km south of Ponce, Puerto Rico. The activity continued with the main sequence starting on December 28 th , 2019 comprising more than 10 Mw>=5.0 earthquakes and the largest, Mw6.4, normal faulting earthquake on January 7 th 2020. The island of Puerto Rico (Fig. 1) is located in a zone of convergence between the North American (NAP) and Caribbean (CP) tectonic plates (Huérfano et al., 1994). North of the 80 island, the NAP is subducting obliquely beneath the CP at the rate of (20 ± 0.4) mm/yr (DeMets et al., 2010) along the east-west striking Puerto Rico trench. Oblique subduction results in strain partitioning and favours microplate tectonics along the eastern Great Antilles island arc (Jansma et al., 2005). The Puerto Rico microplate (PRm) is part of the CP and is moving ENE relative to NAP (Jansma et al., 2000). Rare moderate size earthquakes occurred in and near Puerto Rico in the 20 th century and earlier (Reid et al., 1919, Doser et al., 2005. The majority of instrumental earthquakes (Mw>=3.0) is 95 located north of the island, along the NAP subduction zone. Seismically most active are the SW and SE regions of Puerto Rico (Huérfano et al., 2005). The shallow seismicity of southern Puerto Rico is related to the internal deformation within the Puerto Rico fault zone in the NW-SE direction, while the intermediate depth seismicity occurs above a north dipping zone beneath the island (Huérfano et al., 2005). The tectonics of the southern part of the island is 100 partially controlled by E-W oriented strike-slip faults; South Lajas, Punta Montalva and Salinas faults, which form a more than 80km long fault system. Recent paleoseismologic studies showed that some of the strike-slip faults in the SW of Puerto Rico ruptured in the last 10.000 years and pose a considerable seismic hazard. The South Lajas fault experienced two surface rupturing events, around 5.000 years ago (Prentice and Mann, 2005), while the Salinas fault 105 ruptured twice in the past 10.400 years (Piety et al., 2018). We study the geometry (Fig. 1) of the faults activated in the SW of Puerto Rico during the 2019-2020 Indios sequence by relocating the sequence and performing the kinematic sourcerupture modelling of the five largest events (5.6<Mw<6.4). The relocated events and kinematic rupture inversion reveal a sequence that activated a parallel set of WNW-ESE striking strike-110 slip faults, connected by an orthogonal, NNE-SSW strike-slip fault and oblique, NE-SW striking normal faults. We adopt the matched filter event detection approach to study the temporal and spatial evolution along the faults and in combination with GPS data analysis. We suggest that the sequence was preceded by a transient deformation that most likely loaded the system. 115

Methods
We obtained an earthquake catalog for the Puerto Rico area (latitude 17.56º to 18.73º, longitude -65.50º to -67.80º) from 2019-01-01 to 2020-05-02, M≥2.0 through earthquake.usgs.gov 120 (Benz, 2017) including corresponding PRSN and USGS P and S arrival times, time uncertainties and first-motions. This catalog contains about 7,000 events. We obtained from IRIS-DMC metadata for permanent and temporary stations used for relocation (within 1º of latitude 38.16º, longitude -117.88º; near the Mw6.5 epicenter) and waveforms for coherence analysis, matched-filter and kinematic source analysis. 125 For the geodetic analysis we use the IGS14 continuous GPS (cGPS) position time-series of the Puerto Rico Island stations provided by the repository of the Nevada Geodetic Laboratory (Blewitt et al, 2018).

Precise, absolute earthquake locations 130
We obtained absolute relocations for the study area using source-specific, station travel-time corrections (SSST; Richards-Dinger & Shearer, 2000) and high-precision relocation based on waveform coherency between events. This approach produces enhanced relative location and clustering of events, and allows analysis of model-and data-dependent error in the hypocenters (Lomax, 2020). We performed absolute earthquake location with the NonLinLoc algorithm 135 (Lomax et al., 2000(Lomax et al., , 2014; NLL hereafter), which uses efficient global sampling algorithms to obtain an estimate of the posterior probability density function (PDF) in 3D space for absolute hypocenter location. The location PDF provides a complete description of likely hypocentral locations, includes comprehensive uncertainty information, and allows robust application of waveform coherency relocation. Within NLL, we used the equal differential-time (EDT) 140 likelihood function (Font et al., 2004;A. Lomax, 2005A. Lomax, , 2008Lomax et al., 2014;Zhou, 1994), which is very robust in the presence of outlier data caused by large error in the arrival-times or predicted travel-times. We performed all locations using a search range in depth from 1 to 60 km depth. The shallow limit is a prior constraint required to suppress unreasonably shallow, maximum likelihood depths for numerous events which exhibit double PDF solutions in depth. 145 We used a smoothed, minimum 1D P wave velocity model (Fig. A1) or Puerto Rico and the U.S. Virgin Islands (Huérfano et al., 2010) to seed NLL SSST relocation. We used a finitedifferences, eikonal-equation algorithm (Podvin & Lecomte, 1991) to calculate travel-times for P phases for each station, and obtained S phase travel times from these P times through Vp/Vs=1.71. 150 We relocated events in the Puerto Rico catalog in two stages. First, starting with the initial NLL locations without station corrections, we iteratively generated SSST corrections which vary smoothly within a 3D volume to provide a source-position dependent correction for each station and phase type. We used smoothing distances of 32, 16, 8 and 4 km. Only P and S arrivals with residuals of ≤2.0 sec for relocated events meeting minimum quality criteria are 155 used for update at each iteration. See supplementary material Methods A1 for more details. We relocated the full catalog using the 5 km smoothing-length, SSST corrections (Fig. A2a). In the second relocation stage we greatly reduced absolute location error by combining location information across events based on waveform coherency between the events. This absolute coherency relocation is based on the concept that if the waveforms at a station for two or more 160 events are very similar (have high coherency) up to a given frequency, then the distance separating these "multiplet" events is small relative to the seismic wavelength at that frequency and the events represent stress release on the same, small fault patch (Geller & Mueller, 1980;Poupinet et al., 1982Poupinet et al., , 1984.
The NLL coherence relocation for a target event is a stack over 3D space of the event's SSST 165 location PDF and the SSST PDF's for other events, each weighted by the waveform coherency between the target event and the other event. Unlike differential-time based relative location, absolute coherency relocation requires waveforms from as few as one station, allowing precise relocation for sparse networks, and for foreshocks and early aftershocks of a mainshock sequence or swarm before temporary stations are installed. See the Supplementary material 170 Methods A2 for more details.
The final coherence locations are shown in Fig. A2b. The median formal errors (e.g. with velocity model fixed and assuming only Gaussian, aleatoric error) for the NLL, NLL-SSST and final, NLL-SSST-coherence locations in both models are listed in Table 1 We used waveform data from the 18 closest available stations of the PR seismic network from the beginning of 2019 until the end of January 2020. The waveform data were filtered between 2 and 8 Hz, where the peak energy of most events is expected and their signal to noise ratio is highest, and downsampled to 20 Hz to speed up processing. We prepared a template catalogue of all initial catalog events with semi-major axis <= 5km for 185 the detection of additional, similar earthquakes. The templates were cut from 2.5s before to 2.5s after manually picked S phase arrival times on all available channels and stations. We visually inspected all the templates and removed those of low quality or those that had problems with missing or bad data. To detect new events, the templates were correlated with the continuous waveform data. The 190 sliding-window cross-correlation function (CCF) for each template was calculated with 1 sample step. CCFs were calculated for each station and channel for individual templates and then stacked to form a mean daily CCF trace. For a positive detection we set a threshold of 12*median absolute deviation of the daily trace. After a successful scan of the data we only selected events with inter-event times larger than 3s in order to not include the same event 195 multiple times due to detections from multiple templates. The candidate for detection inside this time window was that with the highest threshold value. The magnitude of the detected event was calculated as the median value of the maximum amplitude ratios for all channels between the template and detected event with a 10-fold increase in amplitude corresponding to a one-unit increase in magnitude (Peng, et al. 2009). 200

Kinematic source inversion
The Indios sequence included five moderate 5.6<Mw<6.4 events (Fig. A3, Table A1; hereinafter EV#1-EV#5). 205 For kinematic source inversion we used near-field strong motion displacement time series derived from digital three-component accelerograms of the PRSN network (Fig. A3). The displacement data were rotated to a defined cartesian coordinate system and filtered with an acausal, fourth-order Butterworth filter in the frequency ranges between 0.03 -0.4 Hz. The selected frequencies depend on the quality of data and the used 1-D velocity model of area. 210 The (usually up to three). Nine parameters are enough to define each slip patch: five to define the ellipse geometry and four for slip, rupture speed, rake, and rise time. 215 The inversion process consists of sampling the trial rupture models by the neighborhood algorithm (Sambridge, 1999a, b) in widely defined ranges for the parameters. Trial synthetic displacement wavefields at each station were computed with Axitra code (Cotton and Coutant, 1997). The synthetic wavefields were compared to observations using the Spudich & Miller (1990) cost function which defines a waveform fit in percent. After sufficient iterations of 220 sampling trial models by the adopted algorithm (usually up to 700), we reached convergence to the final model. More details on the inversion procedure are given in Momeni et al., (2019). We used the hypocenters obtained in section 2.1 as initiation points of ruptures (Table 1). Considering the possible errors in hypocenter location and origin time, we allowed the hypocenters to shift ±1 km along both strike and dip on the fault plane. First, we searched 225 among the nodal planes of the GCMT reported focal mechanism (see Data and Resources) for the geometry of each rupture that provides a rupture model with best waveform-fit to the observations. Selected rupture geometry for all the events are provided in Table 1. After selection of geometry for each event, we ran 10 independent inversions each having slightly different model parameters space. With this we investigate different possible rupture models that can provide the same waveform-fit to the observations. The 10 final rupture models 235 for each event are reported in the Supplementary Material A3 and the preferred rupture models for these events are chosen based on closeness of their scalar seismic moments to the GCMT and USGS point source inversion results (Fig. 2).

GPS 240
We analysed the position time-series (see Data and Resources) following the procedure described in Barzaghi & Borghi, (2018). By considering the temporal correlation among the data and using the least-square estimation method we determined the linear trends (tectonic velocity), seasonal signals, discontinuities due to the station equipment and reference frame 245 changes or seismic events. In Fig. 4a we report the estimated velocity vectors of the cGPS stations for the planimetric components and the velocities with respect to the NAP. Vertical velocities are mostly positive and are reported in Fig A10. The strain-rate principal axes shows a NS stretching especially across the Great Southern Puerto Rico Fault Zone (GSPRFZ, Piety et al. 2018) (Fig. 4c). All of the cGPS stations are characterized by spatial correlation signals, 250 also known as the common mode error (CME  (Table A6). We focused our investigation in the period from 10 th of September 2019 to 27 th December 2019 because it represents a continuous time range with small data gaps and 260 was characterized by increasing seismicity rate. To increase the signal-to-noise ratio of the data and fill the small gaps, the residual time-series were temporally filtered testing two independent methods: the least-square collocation method (LCS) (Borghi et al. 2009;Borghi et al. 2016) using the covariance functions obtained by the analysis of the temporal correlation, and the moving average method (MA) with data windows of two weeks. As described in the 265 supplementary material the two filtering methods provided analogous results (Fig. A11). PCA and BSS methods, like FastICA and variational Bayesian Component Analysis (vbICA) were applied on the residual time-series for transient signals research and will be discussed in the following chapter. 270

Geometry of the activated fault system
The kinematic source inversion of the five largest (5.6<Mw<6.4) earthquakes (Fig. 2, Figs.  275 A4-9) and precisely located fore-and aftershocks show that at least five faults were activated. EV#1 nucleated at the depth of 4.5±1km. The rupture evolved mostly toward up-dip and slightly to the W with an average rupture velocity of 1.8km/s (Vs = 3.1) for ~2.6s. It released a total scalar seismic moment of 5.7E17Nm which is comparable to the GCMT and larger than USGS results (4.5E17Nm and 3.16E17Nm, respectively). The rake is well resolved (-21°) 300 showing a left-lateral strike-slip mechanism. A maximum slip of 1.8m is observable at the depths of 3 to 4km, between 0.7s and 1.7s after the nucleation. EV#2 nucleated at the depth of 8.6km. The rupture evolved on an NE striking (44 °) SE dipping (58°) fault plane with an average speed of 1.55km/s. It reached the surface ~8s after the nucleation (Fig. A6b). The rupture occurred in 10s and released a total scalar seismic moment 305 of 5.6E18Nm, close to the GCMT and USGS results (4.72E+18Nm and 5.04E+17Nm, respectively). It showed a reasonable waveform fit of 75%. Average rake was -124° revealing a normal faulting mechanism with a right-lateral strike-slip component. EV#3 rupture started at the depth of 6.8km. The slip evolved toward shallow depths and slightly to the NE on a NE striking (45°) SE dipping (52°) fault. The whole rupture occurred in 3.3s 310 and released a total scalar seismic moment of 2.74E17Nm. The waveform fit was still reasonable (82%) (Fig. A7c). EV#4 nucleated at a depth of 13km. Its slip evolved mostly down-dip and towards E on an ~E striking, N dipping fault plane with a rupture velocity 2.5km/s, lasting 2.4s with a scalar seismic moment of 6.3E17Nm, larger than the GCMT and USGS results (5.5E17Nm and 3.63E17Nm, 315 respectively). This model has an 84% waveform fit. The rake is 11°, showing almost pure leftlateral strike-slip mechanism. On the 11 th January, EV#5 nucleated at the depth of 19km. Its rupture evolved up-dip and mostly to the N on a ~N-S striking (195°), West dipping (67°) fault plane with an average speed of ~2.8km/s that lasts for ~3.5s. It released a total scalar seismic moment of 1.4E18Nm. This 320 model has a waveform fit of 89%. The obtained rake was -153° showing a right-lateral strikeslip mechanism with a rise time of ~1.6s.

Temporal and spatial evolution
Using well located earthquakes as templates we were able to detect additional events, which gave a better understanding of the sequence development. The sequence began as a low magnitude (M<3.5) swarm, 14km south of Ponce, Puerto Rico in 330 July 2019, ending in January 2020. The newly detected earthquakes in this swarm were temporally split into 3 clusters (Fig. A13) with the most and largest earthquakes in the third cluster.
The main sequence initiated along the Punta Montalva fault (Fig. 3a, b). Prior to the sequence, detected seismicity along the fault was spread over time and we were unable to detect more 335 events than reported in the original catalog. The sequence started on the 28 th December 2019 with a Mw4.7 earthquake, followed a few hours later by a Mw5.0 earthquake and its aftershock sequence. After the Mw5.0 the aftershocks slowly extended both in WNW and ESE direction (Fig. 3a) with the migration velocity towards WNW slightly higher than towards ESE. The WNW side of the Punta 340 Montalva fault was more productive in terms of aftershock than its ESE extent. The ESE migration continued until the 6 th of January 2020 when the largest, Mw5.8 earthquake took place in the ESE portion of the Punta Montalva fault. The Mw5.8 earthquake was preceded by a foreshock sequence of its own in a 2 km wide area of elevated seismicity less than 1 km away from its hypocentre (Fig. A17). In the aftermath of the Punta Montalva fault earthquake, the 345 aftershocks were distributed all along the previously activated fault segment. Few hours after the Punta Montalva fault earthquake an aftershock with Mw4.9 happened at the most ESE part of the Punta Montalva fault. The Mw4.9 earthquake was preceded by an acceleration phase (Fig. A18) at the edge of the activated Punta Montalva fault that stopped with the earthquake. The migration with unchanged velocity continued towards ESE until 13 th of January 2020 when 350 an acceleration phase (Fig. A19) started in the vicinity of a Mw5.2 earthquake that ruptured the most ESE segment of the reactivated part of Punta Montalva fault. On 7 th of January 2020 the Mw6.4 mainshock of the sequence occurred on a normal fault between the Punta Montalva fault and the northern strike-slip fault (Fig. 3b). It was followed by a Mw5.8 normal faulting earthquake NE of its hypocentre. Following the Mw5.8 event, a 355 Mw5.7 strike-slip earthquake ruptured the northern strike-slip fault. The aftershocks along the normal fault were confined towards the north with the area where the fault intersects with the northern strike-slip fault and towards the south, with an intersection with Punta Montalva fault. The aftershock sequence on the normal fault was less energetic than on the northern strike slip fault but this changed on the 8 th January, 2020, when a Mw4.7 earthquake took place on it (Fig.  360 3b) On 11 th January, an orthogonal fault (Fig. 3b) was activated with a Mw6.0 left-lateral strike slip earthquake. Towards the N, aftershocks were confined by the northern strike slip fault, while towards the S, aftershocks extended past the Punta Montalva fault. 365

Geodetic transients
The PCA method applied on the filtered residuals showed that the most important contribution was described by the first two principal components that describe around 80% of the variance. 375 The northern and vertical components are characterized by a discontinuity between October and November 2019, confirmed by Bayesian inference test (Borghi et al., 2012(Borghi et al., , 2016 on the 21 st of October, 2019 (day of the year 294 ± 2 ) with a total probability of 99%. A posterior probability related to the Bayesian estimation of the discontinuity is presented in Figure A12 where we observe the effect of the filtering method used for the detection of discontinuity 380 epoch. PCA and vbICA show very similar behavior in the first principal (PC1) and independent components (IC1) (Fig. A13) and confirm the individuation of a discontinuity in the northern and vertical components. The eastern component does not show a sharp discontinuity (Fig. A4) but rather a positive linear trend that is correlated with the cumulative number of earthquakes in September 2019 (Fig. A14). PC1 could be interpreted as CME signal that should be removed 385 from the time-series, however some considerations suggest that the PCs and ICs reported in Figure A13 represent a geophysical signal present all over Puerto Rico Island in October 2019. Since all the cGPS contribute to the first components (Fig. A13, A15) we suppose that this corresponds to the complex tectonic signal. Moreover, the observed discontinuity temporally falls in the third seismic cluster of Caja swarm, characterized by the most and largest 390 earthquakes (Fig. A14) as already mentioned in Chapter 3.2.
Due to the discontinuity corresponding to the period of increasing seismicity, we analysed this period and computed the displacement field related to this event for each station by least square estimation. The displacement field vectors point towards the south (Fig. 4b) while at the same time, stations undergo uplift (Fig. A16). We also note that the principal axes of the 2d strain 395 tensor are not aligned to the principal axes of the strain-rate tensors due to the tectonic velocity of the stations, but rather represent a divergence behaviour in the NE-SW direction (Fig. 4c)

405
The geometry defined by the relocated events and inversion of the well recorded Mw5.6+ earthquakes on strike-slip and normal faults highlights a very complex system of faulting. We observe that the northern and Punta Montalva faults form a system of parallel strike-slip faults while the normal fault, that hosted the mainshock of the series, forms an oblique structure between them. Additionally, the two parallel strike slip faults are connected by an orthogonal 410 fault deepening from south northward. Its northern edge is bounded by the northern strike-slip fault, while towards the south, it continues past the Punta Montalva fault. These observations show that the geometry of faulting within this sequence is controlled by a larger structural feature -the subduction of CP underneath Puerto Rico as presented in Figure  5.