Shallow Slow Slip Events Identified Offshore the Osa Peninsula in Southern Costa Rica From GNSS Time Series

Using new continuous geodetic time series, we identify five shallow slow slip events (SSEs) offshore and beneath the Osa peninsula in southern Costa Rica. An early event was detected by one station in 2013, and two events occurring in close succession in both 2018 and 2022 were detected by multiple stations, indicating a preliminary recurrence interval of ∼4–5 years. While SSEs have been observed to the northwest at Nicoya, this is their first documentation in southern Costa Rica. Modeled slip distributions of the 2018 and 2022 events indicate they likely ruptured the same or overlapping patches of the plate interface, near the trench, updip of the 1983 Mw 7.4 Osa event. Immediately offshore, estimated cumulative slip from the 2018 and 2022 events is sufficient to close the slip deficit from tectonic loading over the recurrence interval, potentially limiting the magnitude and spatial slip distribution of future large ruptures.

10.1029/2023GL104771 2 of 12 et al., 2014;Kato et al., 2019).This paper focuses on the Osa peninsula in Costa Rica, which along with the Nicoya and Burica peninsulas, extend into the Pacific Ocean nearly to the Middle American Trench and allow for near-trench (often <40 km) terrestrial observations (Figure 1), lending insight into shallow subduction zone processes.
The tectonics of Costa Rica are driven by the northeastward subduction of the Cocos plate and the associated seismicity and volcanism.Near the Osa Peninsula in southern Costa Rica, the Cocos plate subducts to the northeast beneath the Panama microplate at ∼80 mm/yr relative to a stable Caribbean reference frame (DeMets et al., 2010), oriented nearly perpendicular to the Middle American Trench (Figure 1), making it one of the more rapidly converging subduction zones in the world (Schellart et al., 2011).The subduction of the relatively buoyant and thick crust of the Cocos ridge dominates the onshore tectonics, with prior geodetic studies interpreting regional kinematic variations as escape tectonics on either side of the incident Cocos ridge, yielding something akin to forearc sliver motion as well as elevated mechanical coupling at the trench (Kobayashi et al., 2014;LaFemina et al., 2009).
In northern Costa Rica, shallow SSEs have been identified offshore and beneath the Nicoya peninsula (Baba et al., 2021;Brown et al., 2005;Davis et al., 2011Davis et al., , 2015;;Jiang et al., 2012Jiang et al., , 2017;;Outerbridge et al., 2010;Protti et al., 2004Protti et al., , 2005)), near the epicenters of four large (magnitude 7+) earthquakes in 1853, 1900, 1950, and 2012(Protti et al., 2014)).At the Osa peninsula, at least three large (magnitude 7+) earthquakes have occurred in 1904, 1941, and 1983(Adamek et al., 1987;;Tajima & Kikuchi, 1995), with current observations supporting an approximately 40 year recurrence interval.Prior to this study, no SSEs had been documented in this region.The presence of shallow SSEs would therefore shed light on both potential natural hazards in the region and frictional properties of the Middle American Trench.Adamek et al. (1987)."1" is the 20 December 1904 event, "2" is the 5 December 1941 event, and "3" is the 3 April 1983 event.Slab contours are shifted from the Slab2.0model (Hayes et al., 2018) to better align with regional shallow seismicity.Triangles represent regional GNSS stations, with colored stations corresponding to the colored time series plots in Figure 2.

GNSS Data and Processing
To anticipate and better quantify the potential and risks of an impending Osa megathrust earthquake, in recent years the Observatorio Vulcanológico y Sismológico de Costa Rica (OVSICORI-UNA) has been densifying the regional seismic network and establishing a continuously operating Global Navigation Satellite System (GNSS) network in and around the Osa Peninsula.In addition the Insituto Geografico Nacional (IGN) (https://gnss.rnp.go.cr/SBC/) and the Insituto Costarricense de Electricidad (ICE) also operate several continuous GNSS stations in the area.The proximity of these networks to both the offshore trench (Figure 1), and the shallow plate interface (<40 km) (Kyriakopoulos et al., 2015) provide a unique opportunity to use terrestrial methods to study shallow subduction processes including locking and strain accumulation and release, with implications for the entire megathrust.In this study, we use the new GNSS network to identify and describe the first observed shallow SSEs occurring at the Osa Peninsula in southern Costa Rica.
The network of continuously operating GNSS stations was first established in 2010 on and around the Osa Peninsula, with the primary expansion occurring between 2016 and 2018, and additional stations installed in 2022.The majority of these stations are operated by OVSICORI-UNA, two stations (NEIL and SAGE) are operated by the IGN, and others (RC01, PTJI, PN01, PZEN, and BUAI) are operated and maintained by the ICE.We processed the GNSS data using the GAMIT/GLOBK software suite (Herring, King, et al., 2015;Herring, Floyd, et al., 2015), using final orbital products from the International GNSS Service and the Vienna mapping function of Boehm et al. (2006) to correct for tropospheric delays following the processing approach described in Herring et al. (2016).The station position estimates were then adjusted using a Kalman filter and on average 10 stable sites in North America, South America, and the Caribbean to tie the network into the 2014 International Terrestrial Reference Frame (ITRF2014) (Altamimi et al., 2016).Final time series position estimates were then rotated into the stable Caribbean reference frame defined by Altamimi et al. (2012).

SSE Identification and Offset Estimation
We identified SSEs using time series of daily position estimates spanning from 2010 until late 2022 from 28 regional GNSS stations.The processed time series were prepared for examination using a multi-step process.First, we removed outliers, defined as daily measurements with formal errors that are greater than two standard deviations away from the mean position error.Next, we removed non-tidal atmospheric loading using the estimates from the models calculated by Dill and Dobslaw (2013) (Figure S1 in Supporting Information S1).The time series were then cross-referenced with the regional earthquake catalog to check for any offsets related to nearby seismicity and visually inspected for any other transient signals that need to be accounted for.When apparent, instantaneous offsets from earthquakes were removed from the timeseries.We do not account for longer-term post-seismic deformation, primarily because the only event with a measurable offset was an M w 6.1 event on 17 August 2018, and there is no identifiable afterslip decay signal following this event.Finally, we simultaneously estimated inter-SSE deformation as a linear trend, seasonal signals as annual and semi-annual sinusoids, and slow slip offsets as step functions after removal of data points within the period of the transient signal (see Text S1 in Supporting Information S1 for further details).

SSE Modeling
We used the estimated offsets for each SSE to invert for the slip distribution assuming Okada elastic dislocations using GTDef (e.g., Feng et al., 2012;Murekezi et al., 2020).We used the Slab2.0interface geometry (Hayes et al., 2018), which we modified slightly to better align with regional seismicity (Figure S2 in Supporting Information S1).Regions beyond the boundary of Slab2.0 were extrapolated from the established surface.On the downdip edge and sides of the model domain, we constrained slip to be zero (Figure S3 in Supporting Information S1).At the trench we imposed a free surface to allow for slip to the edge of the model domain (Figure S3 in Supporting Information S1), primarily to enable the inversion to better fit observed vertical displacements, though synthetic tests indicate that slip at the trench is best resolved at the far southeast tip of the Osa peninsula (Figure S4 in Supporting Information S1).Additionally, we only inverted for slip using observations from stations that either observed offsets or are located very close to those that observed offsets.Data from relatively far away sites with no clear offset were excluded primarily to prevent deeper slip patches from showing up in the inversion from excess noise.We also constrained the inversion to allow slip to occur with a rake of 60-120°.We then calculated the resolution spread of each model using the method of Kyriakopoulos and Newman (2016) to determine where we could reliably estimate slip distribution.Both horizontal and vertical displacements were used in our inversions, except in one model (late 2022) where some vertical offsets are poorly constrained due to large seasonal signals and thus only horizontal offsets were used.We chose our preferred models and associated smoothing factors based on L-curves of model misfit versus model roughness (Figure S5 in Supporting Information S1) (e.g., Chen et al., 2009;Du et al., 1992;Feng et al., 2012).We then conducted checkerboard resolution tests to visually determine a region where our results are robust (Figure S6 in Supporting Information S1).Lastly, we also constructed several synthetic tests given various input slip distributions to fully characterize where in the model domain is slip resolvable, as well as the implications of different localized slip patches of varying spatial size (Figure S4 in Supporting Information S1).

Results
Time series from more than five GNSS stations exhibit temporally correlated transient signals indicative of four distinct SSEs (Figure 2).An additional transient is apparent at one station in 2013 (Figure S5 in Supporting Information S1), indicating that an earlier SSE is likely to have occurred.

A Potential SSE in 2013
A potential SSE in 2013 was only recorded at one station, RIOS (Figure S7 in Supporting Information S1).The fact that there was only one observation is the result of the GNSS network being relatively sparse in the region at the time.The observed offset, ∼20 mm at RIOS oriented toward the trench (Figure S7 in Supporting Information S1), indicates primarily thrust motions and is similar in magnitude to the offset observed at RIOS in both early 2018 and mid 2022.This indicates that the magnitude of the 2013 SSE and slip distribution are likely similar to the later events, though more detailed modeling of slip distribution is not possible with the available observations.

An Early 2018 SSE
The first event recorded by multiple stations occurred in early 2018 and lasted for roughly one month, with velocity changes recorded at four different GNSS stations.The SSE was observed nearly simultaneously along a trench-perpendicular transect (PIRO, PJIM, and PNPB) with primarily trench-normal displacements that decrease landward (∼60 mm at PIRO, ∼25 mm at PJIM, and ∼2.5 mm at PNPB; Figures 2a and 2c).The SSE was also recorded at RIOS (∼20 mm offset), although the transient began later in time than the other three stations indicating a northwestward propagation of slip.Based on the estimated slip distribution this event occurred ∼20 km updip of a contemporaneous seismic swarm containing ∼80 recorded events during the span of the SSE (Figure 3a) and has an estimated maximum slip of ∼30 cm, equivalent to a M w ∼6.7 earthquake.

A Late 2018 SSE and a M w 6.1 Earthquake
The second event in 2018 is smaller in magnitude than the first and occurred over a span of about one month.Five stations (PNBC: ∼20 mm, PIRO: ∼25 mm, PJIM: ∼15 mm, PTJI: ∼15 mm, and PNPB: ∼5 mm) show an easily identifiable transient in the trench-normal component (Figure 2a), and a muted trench-parallel component (Figure 2c).The displacement appears to have occurred nearly simultaneously with the transient signal observed at PTJI and PJIM occurring a few days prior to the others.Having a maximum estimated slip of ∼15 cm (Figure 3b), this event is equivalent to a M w ∼6.5 earthquake.However, as this SSE occurred approximately 1 week after a M w 6.1 earthquake, the magnitude and spatial distribution of slip may be biased due to the presence of afterslip in some of the offset calculations, particularly in those closer to the epicenter of the M w 6.1 earthquake and its aftershock sequence (Figure 3b).The M w 6.1 earthquake likely occurred within the subducting slab as the hypocenter (∼22 km deep) is below the interface of our modified Slab2.0 geometry (Figure 3b) and the nodal plane closest to the interface is too steeply dipping (∼40° nodal plane vs. ∼20° interface).The aftershock sequence consists of a M w 5.2 event two days after the mainshock as well as 28 events with M w ≥ 3 in the month following the mainshock.

An Early 2022 SSE
Four years later in 2022, we observe another two SSEs.The first occurred in early 2022 and lasted for about one month.This SSE is most apparent in the trench-normal observations at PIRO, PTJI, and PNPB (∼40, ∼10, and ∼10 mm, respectively) with the offset magnitude decreasing away from the trench (Figure 2b).A small trench-parallel (to the northwest) component is also apparent at PIRO (Figure 2d).Modeled slip inversions indicate a maximum estimated slip of ∼15 cm (Figure 3c), broadly equivalent to a M w ∼6.5 event.

A mid 2022 SSE
The final SSE we observed occurred about 1 month later in April to May of 2022.This event was first observed at PLAN and RIOS, two stations located in the northwest portion of the Osa Peninsula, with transients at stations to the southeast (PIRO and PTJI) occurring later in time, indicating southeastward slip propagation (Figure 2b).While all observed transient signals in the geodetic time series have the largest magnitude in the trench-normal component (PLAN: ∼35 mm, RIOS: ∼30 mm, PIRO: ∼25 mm, PTJI: ∼15 mm), there is some apparent slip in the trench-parallel component directed to the northwest (Figure 2d).Modeled slip inversions yield estimates of maximum slip of ∼15 cm in a relatively large region just offshore the Osa Peninsula, equivalent to a M w ∼6.7 event.This event differs from previous events in that the observations indicate a larger area of shallow slip along the trench that includes the region that ruptured in all 3 previous events and a region to its northwest, which is reflected in the larger estimated magnitude.

Synthetic Models and Model Resolution
We generated synthetic models (Figure S4 in Supporting Information S1) to determine the robustness of our spatial slip distribution estimates.The results indicate that while there are some notable gaps in the network coverage, the region we determine as having peak slip in all models is within the limits of detection.In all synthetic models with slip at the trench just offshore the southeastern tip of the Osa Peninsula, input slip distributions near the southeastern tip of the peninsula are reasonably resolved (Figures S4a-S4b, S4h-s4i in Supporting Information S1), further supported by our checkerboard tests (Figure S6 in Supporting Information S1).Those with input slip distributions that continue along the trench to the northwest or smaller localized patches to the northwest (Figures S4e-S4i in Supporting Information S1) show poor recovery of the full input slip distribution.Additionally, we can reasonably resolve the depth of input patches on the interface when it is near the southeastern tip of the Osa peninsula (Figures S4a-S4D in Supporting Information S1).Deeper patches of slip in our inversions (Figures 3a and 3b) are also resolvable in our synthetic models (Figures S4k-S4l in Supporting Information S1), though the small calculated offsets of stations surrounding these slip patches indicates that these may be due to observational noise or other signals (i.e., a small amount of afterslip following the M w 6.1 event) that are difficult to isolate.Lastly, our synthetic models in all cases underestimate the slip magnitude (Figure S4 in Supporting Information S1), indicating that our previously mentioned estimated slip magnitudes are likely also underestimated.However, this is largely expected from elastic models due to a combination of spatial smoothing and filtering through the modeled elastic crust.To summarize, while our estimation of the approximate location of slip in the region immediately off the southeast coast of Osa peninsula is robust, the amplitude of peak slip and spatial extent of slip are more poorly resolved, likely underestimated and overestimated respectively.

Underlying Mechanisms
While it is difficult to determine the specifics of the mechanisms responsible for shallow SSEs in this region, commonalities with other regions lend clues.In many places, such as New Zealand and Japan, SSEs are often attributed to subduction of seamounts that have allowed for the subduction of water rich sediments, as well as the introduction of geometric, lithological, and rheological heterogeneities (Barnes et al., 2020;Bell et al., 2010;Todd et al., 2018;Wang & Bilek, 2014;Yokota et al., 2016).Similarly, SSEs have been observed near the subduction of ridges, namely at Isla de La Plata in Ecuador, near the subduction of the Carnegie ridge (Segovia et al., 2018;Vallée et al., 2013).In southern Costa Rica, while the subduction of the bathymetric high of the Cocos ridge may promote the occurrence of SSEs on its own, the spatial correlation of a large trough, the Cocos ridge central graben (Gardner et al., 2013;Walther, 2003), and the region of maximum slip in all SSEs (Figure 3) indicates excess sediment (compared with the surrounding ridge) is being subducted in this region and that compaction of water rich sediment may allow for low effective normal stresses on the interface, promoting aseismic slip (Kodaira et al., 2004;Peng & Gomberg, 2010;Saffer, 2017).
Another interesting observation is that the SSEs observed in 2018 and 2022 occurred in pairs that follow in close succession and rupture the same or overlapping regions of the plate interface (Figures 2 and 3).In 2018, the two events are separated by ∼5 months, while in 2022 the events are only separated by ∼2 months (Figure 3).While details on the slip distribution and propagation are limited by our model resolution and available data, it is possible that the region of maximum slip did not fully rupture in the first event in either 2018 or 2022, but instead achieved a partial stress drop with the interface primed for further slow slip.Alternatively, the first SSE could have effectively ruptured only one of multiple SSE generating "slow asperities" in the region, effectively loading a nearby asperities that are then more likely to slip in the period following the initial SSE.
While it is tempting to invoke changes in stress conditions from the M w 6.1 earthquake that occurred about one week prior to the apparent onset of the late 2018 event, our stress modeling using Coulomb (Toda et al., 2011) indicates very small stress changes in the region of slow slip (<5 kpa) (Figure S8 in Supporting Information S1).This is significantly less than static stress changes that have triggered SSEs in New Zealand (∼0.1-1.0MPa) (Wallace et al., 2018).However, the lower limit for triggering SSEs is likely small and depends on a number of factors including when in the SSE cycle a stress increase occurs (e.g., Shibazaki et al., 2019), rheological and frictional constraints, existing stress distributions, and pore pressure variations.Additionally, based on the plate interface geometry, the M w 6.1 event is likely an intraplate event within the subducted portion of the Cocos plate.However, we cannot discount triggering from dynamic stress changes (Itaba & Ando, 2011;Wallace et al., 2017), or temporarily elevated pore pressures caused by the earthquake (Khoshmanesh & Shirzaei, 2018;Zhai et al., 2019).In contrast to the slow slip events of 2018, neither of the two events in 2022 are temporally correlated with any sizable seismic event (Figures 3c and 3d), indicating that the observed SSE clustering in time is more likely due to frictional and rheological properties and the pre-existing stress distribution on the shallow megathrust than from seismically induced stress variations.

Common Characteristics
In all observed SSEs, our modeled slip distributions occur in nearly the same or overlapping locations (Figures 4b  and 4c).Only the mid-2022 event breaks this pattern, though only with a slip distribution that occurred over a larger spatial extent.The largest offsets are observed at stations located closest to the trench, indicating shallow slip in the portion of the megathrust updip of the seismogenic zone, with most slip accommodated shallower than 10 km in depth.While deeper aseismic slip patches are apparent in our slip inversion for the early 2018 SSE (Figure 3a) and is theoretically resolvable (Figures S4k-S4l in Supporting Information S1), they are constrained either by one site or small calculated offsets (<5 mm).These could also arise from an actual slip distribution that continues further along strike into a region with limited resolution or noise in our observations, respectively.While presence of slip between the trench and the Osa peninsula is constrained by both horizontal and vertical offsets, uplift is observed at the site closest to the trench (PIRO) in all SSEs, and the largest horizontal offsets are observed at the sites closest to the trench (PIRO and PLAN) (Figure 3), indicating that slip must have occurred in the uppermost portion of the megathrust.However, the lack of additional stations on the southern coast of the Osa peninsula in our synthetic slip inversions indicates that our modeled slip magnitudes are likely underestimated (Figure S4 in Supporting Information S1), and slip distributions may continue further along strike.
While we do not invert for temporal variations in slip, it is apparent that the early 2018 event ruptured from the southeast to the northwest, with a transient being apparent later in time at RIOS (Figure 2a).By contrast, the final event in 2022 has an apparent rupture propagation from the northwest to the southeast.Rough estimates of rupture velocity found using the distance between stations and the onset of the SSE transient indicate a propagation velocity of ∼2 km/day for both the early 2018 event and the late 2022 event, which is slow but within the range of 1-100 km/day for most SSEs (Bürgmann, 2018).Additionally, in 2018, both events were also closely associated with elevated levels of seismicity located downdip from the region of maximum slip (Figures 3a  and 3b), indicating that the region downdip of the SSE is partially coupled.Furthermore, seismicity during the mid 2022 SSE also occurs largely on the downdip margin of the estimated slip patches, supporting the larger spatial extent of our estimated slip distribution (Figure 3d).

Conclusion
We observe three individual slow slip episodes, an episode being defined as any number of slow slip events that are clustered in time and rupture roughly the same patch of the megathrust, here the 2013 event, 2018 events, and 2022 events.These episodes apparently occur every ∼4-5 years, though future observations will be required to more accurately define a recurrence interval.If we are to assume a recurrence interval of 4 years and a Cocos plate velocity of ∼80 mm/yr normal to the trench (DeMets et al., 2010), the amount of slip released in the 2018 and 2022 slow slip episodes is nearly sufficient to close the slip deficit in the region of peak slip at the trench (Figure 4).While this indicates that very little convergence is being accommodated by deformation of the upper and lower plate near the trench, this also implies that a future large earthquake is less likely to rupture the updip region of maximum slow slip, limiting its potential size (Vaca et al., 2018).Additionally, while the apparent release of strain by shallow SSEs may limit the tsunami-genic potential of the megathrust in this region, tsunami risk in Costa Rica is thought to be relatively low (Chacón-Barrantes & Protti, 2011;Protti et al., 2001Protti et al., , 2011)).As SSEs typically rupture conditionally stable regions surrounding seismogenic asperities (Bürgmann, 2018;Xie et al., 2020), the slip distributions of these events may yield 10.1029/2023GL104771 9 of 12 some predictive power in determining where a future earthquake rupture is more or less likely to happen.Furthermore, recent rate and state modeling indicates that SSEs can load seismogenic asperities (Meng & Duan, 2023), and have triggered large ruptures (Dixon et al., 2014;Kato et al., 2019;Uchida et al., 2016), though this is not a requirement (e.g., Voss et al., 2018).As the Osa Peninsula has hosted at least three M ∼ 7-7.5 earthquakes (1904,1942,1983) at nearly a 40 year recurrence interval (Adamek et al., 1987), slow slip events like the ones observed may play a role in loading asperities and/or triggering these and future large earthquake ruptures.(https://gnss.rnp.go.cr/SBC/) and José Ana Gonzalez Carvajal from the Insituto Costarricense de Electricidad (ICE) for their work in collecting and providing GNSS data from stations they operate.We also would like to thank Laura Wallace and an anonymous reviewer for helpful and insightful comments that helped improved the manuscript.Instrumentation and data from OVSICORI-UNA used in this research was supported by Costa Rican Emergency Law 8933 and by Universidad Nacional de Costa Rica through the project 0097-2020 "Sistema de monitoreo geodésico (SiMoGeod) de los volcanes y de la tectónica de Costa Rica."We acknowledge the huge work made by OVSICORI-UNA technical group to install and maintain the GNSS network.This research was supported by the National Research Foundation Singapore under its Singapore NRF Investigator scheme (NRF-NRFI05-2019-0009 to E.M.H.), and by the Earth Observatory of Singapore via its funding from the National Research Foundation Singapore and the Singapore Ministry of Education under the Research Centres of Excellence initiative.This work comprises EOS contribution number 534.

Figure 1 .
Figure1.Overview of the Osa peninsula and surrounding region.Large earthquakes are plotted as numbered stars with locations fromAdamek et al. (1987)."1" is the 20 December 1904 event, "2" is the 5 December 1941 event, and "3" is the 3 April 1983 event.Slab contours are shifted from the Slab2.0model(Hayes et al., 2018) to better align with regional shallow seismicity.Triangles represent regional GNSS stations, with colored stations corresponding to the colored time series plots in Figure2.

Figure 2 .
Figure 2. Horizontal and vertical components of time series of stations that observed the two SSEs in 2018 and two SSEs in 2022.Black dashed vertical lines indicate the start and end of the SSEs from visual inspection while the black solid vertical line indicates the M w 6.1 event.The horizontal components are rotated into trench-normal displacements (panels (a) and (b)) and trench-parallel displacement (panels (c) and (d)).Panels (a), (c), and (e) show time series from 2018.Panels (b), (d), and (f) show time series from 2022.

Figure 3 .
Figure 3. Modeled slip inversions for SSEs in 2018 and 2022.Seismic activity that occurred within the time span of the slow slip event is plotted.Only slip patches in regions that we determine as reasonably well resolved in all SSE's are plotted (black outline in Figure S6 in Supporting Information S1).Black dotted lines are contours of slab depth for the input slab geometry.Observed displacements are plotted in cyan while modeled displacements are plotted in magenta.Panel (a) is the early 2018 SSE, panel (b) is the late 2018 SSE and includes the focal mechanism for the M w 6.1 event that occurred ∼1 week prior to the SSE, panel (c) is the early 2022 SSE, and panel (d) is the mid 2022 SSE.

Figure 4 .
Figure 4.Estimated slip at the trench shows that maximum slip generally occurred in the same overlapping area for each event.Panel (a) shows the GNSS stations used in each inversion for each SSE and their distance to the trench, and indicates that the greatest slip was estimated where we have stations closest to the trench.Panel (b) shows slip at the trench for each modeled SSE.Panel (c) shows the cumulative slip for 2018 and 2022 in comparison with the amount of slip that would be required over a 4 year recurrence interval if the interface was slipping at the plate rate.