Tectonics and Sliver Motion in Costa Rica: Strain Partitioning Constrained by GNSS Velocities Data – Boymond – 2026 – Geochemistry, Geophysics, Geosystems

3.1 The GNSS Velocities

The Observatorio Vulcanológico y Sismológico de Costa Rica (OVSICORI), Instituto Geográfico Nacional (IGN), the Instituto Costariciense de Electricidad (ICE) and UNAVCO operate a network of 154 permanent GNSS stations (Figure 2) to detect ground deformation linked to volcanic and tectonic processes. In accordance with the OVSICORI’s mission of monitoring volcanic and tectonic risk in the country, the stations are distributed all throughout Costa Rica with higher spatial density in the at-risk volcanic areas. Most stations have been acquiring position data since 2015, with some going as far back as 2000. For this study, we use data starting in 2015, marking the end of the post-seismic deformation following the 2012 Mw 7.6 Nicoya earthquake through 2024 (Hobbs et al., 2019) as observed in time series data for GNSS station GRZA (see Figure S1 in Supporting Information S1 for time series, Figure 2 for GRZA location).

The equipment used to ensure accurate, quasi-continuous measurements varies from station to station. In general, a station consists of a GNSS antenna (Trimble Zephyr Geodetic 1, 2 or 3 as well as Tallysman 6500) mounted on a concrete base or structure, a GNSS receiver (Trimble NetR9, Trimble Alloy or Septentrio Polarx5), powered by gel or lithium batteries which are recharged during the day by solar panels.

The GNSS data were processed using the GAMIT/GlOBK software (Herring et al., 2018). Tropospheric and ionospheric effects were corrected using the Vienna Mapping Function (VMF1; Boehm et al., 2006) and final ionospheric maps from the International GNSS Service (IGS). We also used satellite final orbits and relevant earth orientation parameters provided by the IGS. In addition, ocean tide loading corrections were applied using the FES2004 tidal model (Lyard et al., 2006). The time series were manually screened for outliers and corrected from the displacements generated by the largest seismic events and equipment changes using standard equations (Bock & Melgar, 2016). To ensure reproducibility, all equipment-related discontinuities and coseismic offsets were estimated using standard GLOBK protocols and are documented in the Tables S1 and S2 of Supporting Information S1. Then, a Kalman filter was applied to adjust station positions relative to a reference frame consisting of at least 18 stable sites (see Table S3 in Supporting Information S1) and the 2014 International Terrestrial Reference Frame (ITRF2014; Altamimi et al., 2016). The final position estimates from the time series were rotated into the stable Caribbean reference frame CARIB08 (Altamimi et al., 2012). We provide a map of the GNSS velocities of this study in the ITRF2014 reference frame, with velocities from previous studies also converted to the ITRF2014 in Figure S2 in Supporting Information S1. Consistent with our objective to characterize long-term regional deformation rather than absolute coupling at the subduction interface (e.g., Perry et al., 2025), individual Slow Slip Events (SSEs) were not removed as their signals are effectively averaged into the linear rates over the 9-year observation period. We obtain the eastward and northward velocities for each station, as well as their error ellipses (Tables S4 and S5 in Supporting Information S1 for the velocities in ITRF14 and in the stable CARIB08 reference frame). As the study of vertical ground motion and deformation is outside the scope of this study, the vertical velocities will not be used but are presented in Tables S4 and S5 of Supporting Information S1.

3.2 Strain Rate and Velocity Field

The GNSS network provides extensive spatial coverage across all of Costa Rica, spanning 51,000 km². However, the distribution of stations is not uniform. Some regions, such as urban areas and volcanic zones, have a higher concentration of stations, while others, such as the Cordillera de Talamanca and the Northern Limón Basin (NLB), have relatively few stations. This uneven distribution creates a limitation in our understanding of the tectonics, deformation, and geodynamics of these regions. To better constrain deformation over the full extent of Costa Rica, we used an interpolation method adapted from the VISR strain rate interpolation proposed by Shen et al. (2015) to compute the velocities and strain rates in between the GNSS stations, using the GNSS velocities fixed in the Caribbean reference frame (CARIB08—Altamimi et al., 2012). It first evaluates the GNSS stations that are statistical outliers based on a criterion on their studentized residuals. After excluding the identified outliers, the least-squares inversion is then recomputed using weights derived from the GNSS velocity uncertainties, the local density of the GNSS network, and a smoothing distance evaluated for each pixel on the interpolation grid (see Text S1 in Supporting Information S1 for the weighting formulation). The interpolation, the choice of smoothing parameter and the process of evaluating the outliers are all described in more detail in Text S1, Figures S3 and S4 of Supporting Information S1. The study area is a 356-wide by 440-long km area encompassing Costa Rica, from the Nicaraguan frontier in the north to the Burica Peninsula in the south. This zone is discretized in M = 498,057 cells measuring 5 by 5 km. For each cell, we interpolate the local north-south (Figure 3a) and east-west (Figure 3b) velocities as well as its strain rate tensor. From the strain-rate tensor, we computed the two principal strains (Figures 3c and 3d). We also computed the first invariant of the strain-rate tensor (Figure 3c) and the second invariant of the deviatoric strain-rate tensor (Figure 3d). The first invariant is commonly used as a proxy for compression (negative) or dilation (positive). The second invariant quantifies the amplitude of non-volumetric (distortional) deformation; in a 2-D horizontal deformation framework, this component primarily corresponds to shear (Gerya, 2019). We therefore interpret the deviatoric second invariant as a proxy for shear strain (see Text S1 in Supporting Information S1 for invariant and principal-strain computations). During interpolation, local uncertainty is also calculated for velocities and strain rates (see Figure S5 in Supporting Information S1 for uncertainty maps and thresholds applied for visualization), which reflects the precision of the results. This uncertainty is represented by the level of transparency in the figures of this study.

The interpolated north-south and east-west velocity fields were projected onto two sets of key azimuths to facilitate interpretation. First, we projected the velocities parallel (N22°E) and perpendicular (N112°E) to the Caribbean Plate convergence direction (Figures 4a and 4b resp.). In central and southern Costa Rica, the N112°E direction also coincides with the azimuth of the subduction trench. Second, we projected the velocities perpendicular (N45°E) and parallel (N135°E) to the strike of the MAT in northern Costa Rica (Figures 4c and 4d resp.). The N135°E azimuth also corresponds to the strike of the Haciendas-Chiripa Fault System (HCFS). For the remainder of the paper, following standard geophysical conventions, velocities projected along a given azimuth are defined as positive in the direction of that azimuth. When comparing velocities with opposite signs, the relative amplitude refers to their absolute values.

Although the strain-rate field is derived using the VISR methodology (Shen et al., 2015), we note that alternative interpolation strategies (e.g., triangulation-based, basis-function, elasticity-based, and Bayesian methods) can yield different amplitudes and spatial patterns, particularly in regions of sparse data or low strain rates (Maurer & Materna, 2023; Métois et al., 2025). Recent intercomparisons indicate that while VISR may exhibit some amplitude attenuation and increased dispersion at short wavelengths, results are generally robust at the regional scale considered here; nevertheless, small localized features should be interpreted with caution as they may partly reflect methodological artifacts.

4.1 Convergence-Parallel and Perpendicular Velocities

The N-S and E-W velocity components of the GNSS stations were reprojected into a reference frame defined by axes parallel (N22°E) and perpendicular (N112°E) to the local Caribbean subduction convergence direction (N22°E) (vectors in Figures 4a and 4b, respectively). Similarly, the interpolated N-S and E-W velocity fields were reprojected onto the same N22°E-N112°E basis (colormaps in Figures 4a and 4b, respectively).

Overall, GNSS velocities and the interpolated velocity field, both projected along N22°E, show a consistent spatial gradient: velocities are highest near the subduction trench and progressively decrease toward the volcanic arc and backarc. The largest values are observed at the southeastern stations of the Osa and Burica Peninsulas, with PJIM reaching $$33.3\pm 0.3\,\text{mm}.\mathrm{y}{\mathrm{r}}^{-1}$$ and PLAN $$32.3\pm 0.1\,\text{mm}.\mathrm{y}{\mathrm{r}}^{-1}$$ (see Figure 4a for station locations). Stations in the central section and at the southernmost tip of the Nicoya Peninsula display the second-highest velocities, with IND1 reaching $$28.5\pm 0.1\,\text{mm}.\mathrm{y}{\mathrm{r}}^{-1}$$. In the interpolated field, convergence-parallel velocities reach up to $$33\pm 0.1\ \text{mm}.\mathrm{y}{\mathrm{r}}^{-1}$$ on the Osa Peninsula and $$26\pm 0.1\,\text{mm}.\mathrm{y}{\mathrm{r}}^{-1}$$ over Nicoya. In contrast, the lowest velocities occur in the North Limón Basin backarc, where GNSS stations such as COVE ($$6.8\pm 0.1\,\text{mm}.\mathrm{y}{\mathrm{r}}^{-1}$$) and MUEL ($$4.5\pm 0.12\,\text{mm}.\mathrm{y}{\mathrm{r}}^{-1}$$) correspond to interpolated values as low as $$1.3\pm 0.14\,\text{mm}.\mathrm{y}{\mathrm{r}}^{-1}$$. Similarly, the blue curve in Figure 5a further shows that, for a given distance from the subduction trench, interpolated velocities in northwestern and southeastern Costa Rica are systematically higher than in central Costa Rica, where maximum values reach only around $$15\,\text{mm}.\mathrm{y}{\mathrm{r}}^{-1}$$ near the Pacific coastline. The velocity profile shown in Figure 5c indicates that, in central Costa Rica, Caribbean Plate convergence-parallel velocities remain nearly constant across the forearc and begin to decrease only 30 km from the volcanic arc inland. In contrast, the profiles in Figure 5b (southeastern Costa Rica) and Figure 5d (northwestern Costa Rica) show an almost linear decrease in convergence-parallel velocities with distance from the trench. For southeastern Costa Rica, a linear regression of the velocity profile was computed and is shown as the gray dotted line in Figure 5b. Deviations from this linear trend reach a maximum of $$0.6\pm 0.1\,\text{mm}.\mathrm{y}{\mathrm{r}}^{-1}$$ and a minimum of $${-}0.8\pm 0.1\,\text{mm}.\mathrm{y}{\mathrm{r}}^{-1}$$ over a zone comprising the Cordillera de Talamanca and Fila Costeña, corresponding to a total peak-to-peak dispersion of around $$1.5\pm 0.1\,\text{mm}.\mathrm{y}{\mathrm{r}}^{-1}$$.

Arrays of Figure 4b display the component of the GNSS velocities along the N112°E azimuth, which corresponds to the direction perpendicular to the local direction of convergence for the Caribbean plate (N22E°). Incidentally, the N112°E also corresponds to the strike of the MAT in southeastern and central Costa Rica, up to Nicoya, where it bends toward the N135°E azimuth. For this study, negative N112°E-parallel velocities indicate a northwestward displacement, while positive N112°E-parallel velocities indicate a southeastward motion.

The velocities of both the stations and the interpolated field are the lowest along the entire backarc region. The N112°E interpolated field reaching down to $${-}0.4\,\text{mm}.\mathrm{y}{\mathrm{r}}^{-1}$$ in the SLB ($${-}0.9\pm 0.9\,\text{mm}.\mathrm{y}{\mathrm{r}}^{-1}$$ for the BRBR station of the SLB—stations locations in Figure 4b) and $${-}3\pm 0.9\,\text{mm}.\mathrm{y}{\mathrm{r}}^{-1}$$ in the NLB ($${-}0.1\pm 2.2\,\text{mm}.\mathrm{y}{\mathrm{r}}^{-1}$$ for the CAPO station in the NLB). In the southeastern part of Costa Rica, the N112°E GNSS velocities and velocity field are low too, with velocities between $$0$$ and $${-}2\,\text{mm}.\mathrm{y}{\mathrm{r}}^{-1}$$ the velocity field, reaching down to $$0\pm 0.2\,\text{mm}.\mathrm{y}{\mathrm{r}}^{-1}$$ in CRUC on the Fila Costena, and $${-}1.6\pm 0.1\,\text{mm}.\mathrm{y}{\mathrm{r}}^{-1}$$ in PLAN on the Osa Peninsula. As illustrated in the profile of Figure 6, the N112°E interpolated velocities progressively increase northwestward along the arc and forearc, ranging from $${-}3.6\,\text{to}-6.8\,\text{mm}.\mathrm{y}{\mathrm{r}}^{-1}$$ in Central Costa Rica, reaching an average of $${-}7.8\,\text{mm}.\mathrm{y}{\mathrm{r}}^{-1}$$ in the Nicoya Peninsula, and $${-}11$$ to $${-}12\,\text{mm}.\mathrm{y}{\mathrm{r}}^{-1}$$ in the Northwesternmost part of Costa Rica ($${-}11.6\pm 0.3\,\text{mm}.\mathrm{y}{\mathrm{r}}^{-1}$$ in LACR).

4.2 HCFS-Parallel and Perpendicular Velocities

The GNSS station velocities were projected perpendicular (N45°E) and parallel (N135°E) to the strike of the HCFS (N135°E) and are shown as vectors in Figures 4c and 4d, respectively. The N135°E axis also corresponds to the azimuth of the trench in northwestern Costa Rica. Figures 4c and 4d additionally display the interpolated N45°E and N135°E velocity components as colormaps.

At the scale of Costa Rica, the N45°E velocity field and GNSS vectors closely resemble the convergence-parallel (N22°E) GNSS velocities and interpolated field. Velocities are highest near the Middle America Trench (MAT) and decrease progressively with distance from the trench, reaching minimum values in the backarc, down to $$1.8\pm 0.2\,\text{mm}.\mathrm{y}{\mathrm{r}}^{-1}$$ in the NLB. Two zones of elevated velocities are observed on the Nicoya Peninsula, with a maximum of $$21.3\pm 0.1\,\text{mm}.\mathrm{y}{\mathrm{r}}^{-1}$$, and in southeastern Costa Rica, where velocities reach up to $$30.7\pm 0.1\,\text{mm}.\mathrm{y}{\mathrm{r}}^{-1}$$. In contrast, central Costa Rica is characterized by comparatively lower but spatially stable velocities, averaging around $$10\,\text{mm}.\mathrm{y}{\mathrm{r}}^{-1}$$.

The N135°E velocity field and GNSS velocities show a decrease in velocity from the MAT, with low velocities in the backarc reaching down to $$0.1\pm 0.2\,\text{mm}.\mathrm{y}{\mathrm{r}}^{-1}$$ in the NLB, and high velocities in the forearc culminating in the Nicoya (maximum at $${-}16.7\pm 0.1\,\text{mm}.\mathrm{y}{\mathrm{r}}^{-1}$$—negative velocities display a northeastward motion) and Osa Peninsulas ($${-}14.5\pm 0.07\,\text{mm}.\mathrm{y}{\mathrm{r}}^{-1}$$). N135°E velocities are stable in Central Costa Rica around $${-}10$$ to $${-}12\,\text{mm}.\mathrm{y}{\mathrm{r}}^{-1}$$.

The velocity profile of Figure 6 illustrates that, for northeast Costa Rica, the high magnitudes for the N135°E component of the velocity field remain stable all throughout the forearc and part of the volcanic arc but decrease progressively throughout the arc and backarc. The N135°E component of the GNSS station velocities in northeastern Costa Rica is also plotted along the profile (green dots) and displays a characteristic velocity gradient across the fault. We fitted an arctangent function commonly used to describe interseismic deformation across strike-slip faults (Savage & Burford, 1973). The best-fitting model (gray line in profile of Figure 6) indicates a slip rate of $$10.8\pm 0.1\text{mm}.\mathrm{y}{\mathrm{r}}^{-1}$$ accommodated along the fault system and a locking depth of $$7\pm 2\text{km}$$. The location of the fitted fault is consistent with that of the HCFS but is located 6 km to the northeast. The depth distribution of seismic events associated with the HCFS (histogram in Figure 6c) indicates that most earthquakes occur between 10 and 15 km depth.

4.3 Strain Rates

The first invariant of the strain rate tensor (j1) is often used as a proxy value for compression, if negative, or dilation, if positive (Maurer & Materna, 2023), whereas the second invariant of the deviatoric of the strain rate tensor is used as a proxy value for shear (Gerya, 2019) (see Text S1 in Supporting Information S1 for the invariant calculations). We observe in Figures 3c and 3d that the invariants are not homogeneous throughout Costa Rica.

Overall, the first invariant of the strain rate tensor is negative all over Costa Rica, indicating a compressional tectonic regime. Northwestern and southeastern Costa Rica correspond to zones with higher compression, with the magnitude of the first invariant reaching a maximum respectively in Nicoya Peninsula ($${-}20.6\pm 0.3\,\text{mstrain}.\mathrm{y}{\mathrm{r}}^{-1}$$) and on the Talamanca Cordillera ($${-}27.4\pm 0.6\,\text{mstrain}.\mathrm{y}{\mathrm{r}}^{-1}$$). In Central Costa Rica, the first invariant of the strain rate tensor is closer to zero, marking a 90-km wide North South corridor of low compression between northwest and southeast Costa Rica. Here, the first invariant for the strain rate tensor even reaches positive values in a small area centered on the Central Valley ($$4.39\pm 0.5\,\text{mstrain}.\mathrm{y}{\mathrm{r}}^{-1}$$), suggesting a local extensional regime, or at least, as suggested by the principal strain components, N-S to NW-SE compression compensated by E-W to NW-SE extension. The red profile in Figure 5a displays the evolution of the first invariant of the strain rate tensor (j1) along a trench-parallel profile transecting Costa Rica from NW to SE, along the Nicoya Peninsula, Central Costa Rica and the Talamanca Cordillera in the south. Once again, it illustrates that the higher compressive tectonic regime is localized mainly in the northeast and southwest zones. The Central zone, where the CCRDB passes through, corresponds to an area where the first invariant is closer to null. Figure 5b, a SSW to NNE profile transecting the forearc, arc and backarc in southeastern Costa Rica, illustrates that the first invariant of the strain rate tensor remains stable around $${-}14\,\text{mstrain}.\mathrm{y}{\mathrm{r}}^{-1}$$ throughout the forearc and backarc, but peaks to $${-}28\pm 0.7\,\text{mstrain}.\mathrm{y}{\mathrm{r}}^{-1}$$ on the Talamanca Cordillera. This indicates higher compressional strain that correlates to where the peak deviation of N22°E component of the interpolated velocity surface to the linear fit ($$1.5\,\text{mm}.\mathrm{y}{\mathrm{r}}^{-1}$$) is located (blue curves Figure 5b).

The second invariant of the strain rate tensor (j2—Figure 3d) reaches its highest values in the southeast, along the Burica and Osa peninsulas and the Talamanca Cordillera ($$11.6\pm 0.3\,\text{mstrain}.\mathrm{y}{\mathrm{r}}^{-1}$$), as well as in the north along the Northern Limón Basin (NLB), northeast of the HCFS (up to $$11.3\pm 0.2\,\text{mstrain}.\mathrm{y}{\mathrm{r}}^{-1}$$). These high values indicate zones of enhanced shear deformation. Figure 6c shows the variation of the second invariant of the deviatoric strain-rate tensor along a N45°E-oriented profile perpendicular to the HCFS. The second invariant is systematically higher on the northeastern side of the HCFS, highlighting significant shearing within the NLB.

In Figures 3c and 3d, both the first invariant of the strain-rate tensor and the second invariant of the deviatoric strain-rate tensor show that the interpolated fields extend offshore, west and south of the Nicoya Peninsula, where no GNSS observations are available. These interpolations suggest high shear west and south of the Nicoya Peninsula, pronounced compression south of the peninsula, and extension to the west. However, because these features lie outside the region constrained by GNSS data, we do not interpret them further.

The principal strain rate components (displayed as vectors in Figures 3c and 3d) indicate NE-SW to ENE-WSW compression in the forearc and volcanic arc in northeastern Costa Rica. The large contrast between the principal strain-rate components in this region’s backarc results in high deviatoric strain rates, indicative of strong shear deformation. In southeastern Costa Rica, the principal strain rate components are indicative of compression along NE-SW axis. In Central Costa Rica, along the CCRDB, principal strain components indicate compression along the N-S axis and extension along the E-W axis in the forearc and arc. In the backarc, the axis of rotation and extension shifts as principal strain rate components indicate compression along the NE-SW and extension along the NW-SE axis.

5.1 Compression in Southeastern Costa Rica

The near absence of obliquity between the Middle America Trench (MAT) and the Caribbean plate convergence direction inhibits slip partitioning at the plate interface, resulting in a negligible trench-parallel velocity component (Figure 4b). Consequently, the compressional strain indicated by the large negative values of the first invariant of the strain-rate tensor (Figure 3c) is accommodated almost entirely by trench-normal motion (Figure 4a). Trench-normal velocities decrease approximately linearly from the MAT toward the backarc. The largest departure from this linear trend is observed across a region comprising the Fila Costeña and Cordillera de Talamanca. Given the small amplitude of this deviation and the smoothing inherent to the interpolation procedure, this value is not interpreted as a direct measure of the deformation accommodated across these structures and may partly reflect residual interpolation effects. Nevertheless, this region is characterized by a pronounced and statistically significant peak in compressional strain rate, indicating that deformation is indeed focused across the Fila Costeña-Cordillera de Talamanca system. This deformation is interpreted to be accommodated through a combination of reverse faulting and distributed plastic deformation (Alfaro et al., 2023; Sitchler et al., 2007), in a region marked by significant crustal thickening and a steepening of the subduction interface beneath the Cordillera de Talamanca, as inferred from seismic and structural studies and interpreted to reflect the absence of subducting Cocos Ridge buoyancy at this latitude (Dzierma et al., 2010, 2011; Norabuena et al., 2004).

Outside this region, the smooth $$30\,\text{mm}.\mathrm{y}{\mathrm{r}}^{-1}$$ velocity gradient may reflect predominantly interseismic strain accumulation at the plate interface, potentially associated with elastic loading of the forearc (Kobayashi et al., 2014; LaFemina et al., 2009) prior to large Mw > 7 megathrust earthquakes such as those of 1904, 1941, and 1983 documented offshore the Osa Peninsula (Perry et al., 2025). The absence of a marked velocity discontinuity suggests that interplate coupling between the Cocos and Caribbean plates extends beyond the Osa Peninsula, both laterally and downdip, rather than being confined to a narrow, highly coupled patch on and offshore Osa. While the Osa Peninsula may host locally very high coupling coefficients (approaching 100%, Perry et al., 2025), the GNSS velocity field indicates that elastic strain is distributed over a broad region of the overriding plate. This pattern may result either from laterally and downdip extensive coupling on the megathrust, or from a mechanically strong overriding plate that distributes elastic deformation over a wide area rather than localizing it near individual structures. Compared to the lower values observed in central Costa Rica, the compressional strain rates are consistent with enhanced shortening of the overriding plate driven by the buoyant subduction of the Cocos Ridge; in turn, this might drive the northeastward motion of the Panama Block (Kobayashi et al., 2014; LaFemina et al., 2009; Van Benthem & Govers, 2010). However, our results do not show evidence of a trench-parallel push associated with this indentation.

Because our analysis relies exclusively on horizontal GNSS velocities, vertical deformation associated with crustal thickening, uplift, or flexural loading is not captured; incorporating vertical GNSS velocities would provide an important complementary constraint on the partitioning between permanent deformation and interseismic elastic strain accumulation in this region.

5.2 Sliver Tectonics in Costa Rica

This study reveals a clear northwestward increase in the trench-parallel component of velocities across Costa Rica (Figure 5a), from negligible values in the southeastern forearc (Osa and Burica peninsulas) to significantly higher values in the Nicoya Peninsula and northwestern Costa Rica. This first-order pattern, observed consistently across both the forearc and volcanic arc domains, indicates a progressive reorganization of deformation. As the Middle America Trench (MAT) rotates from a N112°E orientation in southeastern and central Costa Rica to ∼N135°E in the northwest, trench-parallel velocities broadly follow this geometric trend but become markedly enhanced in northwestern Costa Rica. This behavior is consistent with increasing strain partitioning of oblique convergence toward the northwest.

Our results refine the along-arc evolution of forearc sliver motion and are broadly consistent with previous studies, which report moderate trench-parallel velocities in the Nicoya Peninsula (Feng et al., 2012; Hobbs et al., 2019; Iinuma, 2004; LaFemina et al., 2009; Lundgren et al., 1999; Norabuena et al., 2004) increasing farther northwest along the Central American margin (Álvarez-Gómez et al., 2019; Kobayashi et al., 2014). However, these comparisons are not intended as point-by-point equivalences, as the reported velocities are derived from data sets spanning different time periods, spatial extents, and processing strategies, and expressed in distinct (though Caribbean-fixed) reference frames.

Similar trench-parallel forearc motion attributed to sliver tectonics has been documented elsewhere along the Middle America Trench (e.g., El Salvador—Álvarez-Gómez et al., 2019; Dewey et al., 2004; Ellis et al., 2019, southern Guatemala—Ellis et al., 2019; Garnier et al., 2021) and in South America (Ecuador and northern Peru—Alvarado et al., 2011; Tamay et al., 2021), with GNSS studies showing that convergence obliquity and interplate coupling can play a major role in controlling such translation (Chlieh et al., 2014; Nocquet et al., 2014; Villegas-Lanza et al., 2016). In this context, the near-zero trench-parallel velocities observed in southeastern Costa Rica provide a key constraint on the driving mechanism. Specifically, the observed northwestward increase strongly favors a drag-driven model (Álvarez-Gómez et al., 2019), in which increasing obliquity between plate convergence and trench orientation promotes partitioning into trench-parallel and trench-normal components. In contrast, a push model alone (Kobayashi et al., 2014), associated with lateral indentation by the Cocos Ridge, would tend to predict a different along-strike pattern of trench-parallel velocities than observed here, which appears less consistent with the observed northwestward increase in trench-parallel velocities. Although a combined push—drag mechanism cannot be excluded (Álvarez-Gómez et al., 2019; LaFemina et al., 2009), our results indicate that strain partitioning driven by increasing convergence obliquity is the dominant control on forearc motion in Costa Rica and farther northwest.

Along obliquely convergent margins, forearc segmentation expressed by alternating embayments and peninsulas has been shown to reflect spatial variations in slab geometry, sediment input, and interplate friction, which modulate strain partitioning and trench-parallel motion (Cubas et al., 2022; Saillard et al., 2017). The Nicoya Peninsula constitutes a major morphological segment of the Costa Rican margin and coincides with enhanced trench-parallel velocities, consistent with this framework.

5.3 The HCFS

In this study, the HCFS is simplified as a single continuous fault because the GNSS network coverage is not locally dense enough to resolve individual faults and fault segments composing the HCFS. Our results show that this amalgamated fault accommodates a $$10.8\pm 0.1\,\text{mm}.\mathrm{y}{\mathrm{r}}^{-1}$$ velocity with a $$7\pm 2\,\text{km}$$ locking depth. The zone of relatively higher second invariant of the deviatoric strain rate tensor located on and northeast of the fault system provides evidence for a shearing tectonic context, indicating that this $$10.8\pm 0.1\,\text{mm}.\mathrm{y}{\mathrm{r}}^{-1}$$ slip velocity can be accommodated by dextral strike-slip motion along the HCFS. This kinematic pattern aligns with the interpretation of other authors of the HCFS as the sliver block boundary in northwestern Costa Rica (Araya & Biggs, 2020; Montero et al., 2017), extending sliver tectonics beyond the forearc to include the volcanic arc and parts of the back arc.

We show that trench-parallel sliver velocities in northwestern Costa Rica range between $$12$$ and $$16\,\text{mm}.\mathrm{y}{\mathrm{r}}^{-1}$$, implying that the HCFS accommodates approximately 63%–83% of the northwestward motion of the sliver block. However, GNSS and seismic network coverage within the NLB is insufficient to resolve how the $$10.8\pm 0.1\,\text{mm}.\mathrm{y}{\mathrm{r}}^{-1}$$ HCFS is partitioned between aseismic creep, seismic slip, and elastic strain accumulation associated with future seismic release. Consequently, we are unable to assess the spatial pattern of interseismic coupling along the HCFS.

The interpolated principal components of the strain-rate tensor indicate E-W compression in the vicinity of the HCFS (Figure 6b), consistent with transpressional deformation along faults striking approximately N135°E. This interpretation is further supported by the focal mechanisms associated with the HCFS (Figure 6b), which display a combination of strike-slip and compressional faulting and the positive flower structures proposed by Araya and Biggs (2020) for the Caño Negro and Upala faults of the HCFS (Figure 1), which are features commonly associated with sliver block transport in transpressional tectonic settings (e.g.,Bénâtre, 2023; Chizzini et al., 2022). Taken together, these observations support the interpretation of the HCFS as a transpressional system accommodating sliver motion. Therefore, evaluating whether this deformation is restricted to the crust or extends into the lithosphere becomes essential. The shallow depth distribution of seismic events associated with the HCFS (histogram in Figure 6c) indicates that most earthquakes occur within the upper crust. Although the maximum depths exceed the locking depth inferred in this study, this discrepancy may reflect the relatively sparse seismic network coverage in the NLB compared to the rest of Costa Rica, which can result in incomplete detection of seismicity and limited precision in hypocentral depth estimates, thereby precluding a robust comparison. Nevertheless, the observed seismicity remains confined to shallow crustal levels, indicating that deformation along the HCFS is primarily accommodated within the upper crust. Taken together, these observations support the interpretation of the HCFS as a crustal-scale kinematic discontinuity accommodating sliver motion, rather than a lithosphere-scale plate boundary. To conclusively resolve the nature of the HCFS, particularly its depth extent, and to discriminate between a true flower structure and a system of closely spaced, mechanically decoupled transpressive faults, seismic imaging capable of resolving deeper structural connectivity is required.

5.4 Central Costa Rica and the CCRDB

In Central Costa Rica, the first invariant of the strain-rate tensor (Figure 3c) is comparatively lower than in northwestern and southeastern Costa Rica and close to $$11\,\text{mstrain}.\mathrm{y}{\mathrm{r}}^{-1}$$. This indicates deformation with little net change in area (hereafter referred to as isochoric deformation), meaning that, in the principal strain-rate frame, extension is balanced by shortening, consistent with predominantly distortional (shear) deformation, in contrast to the compressive deformation observed in northeastern and southwestern Costa Rica. This low-invariant domain forms a rough N-S trending band that separates the northwestern and southeastern regions and coincides with the CCRDB.

Analysis of the principal strain-rate components (Figures 3c and 3d) further reveals a clear contrast between the forearc and backarc domains within this central zone. In both domains, the first $$\left({\dot{\varepsilon }}_{1}\right)$$ and second $$\left({\dot{\varepsilon }}_{2}\right)$$ principal strain-rate components indicate compressional strain $$\left({\dot{\varepsilon }}_{1}< 0\right)$$ associated with extension $$\left({\dot{\varepsilon }}_{2} > 0\right)$$. However, in the forearc, the magnitudes of the principal strain components are smaller than in the backarc. There, $${\dot{\varepsilon }}_{1}$$ is oriented N-S and $${\dot{\varepsilon }}_{2}$$ E-W, whereas in the backarc these axes rotate to ESE-WNW for $${\dot{\varepsilon }}_{1}$$ and NNE-SSW for $${\dot{\varepsilon }}_{2}$$. This rotation corresponds to the geometric bend of the CCRDB.

This pattern indicates elevated deviatoric strain, as confirmed by higher values of the second invariant of the strain-rate tensor in the back arc (Figure 3d), and reflects isochoric deformation. In the backarc, shortening along a WNW-ESE axis is compensated by extension along an NNE-SSW, and in the forearc, N-S shortening is compensated by E-W extension.

We thus propose that the CCRDB acts as a buffer zone between northwestern and southeastern Costa Rica, accommodating a northwestward migration of northwestern Costa Rica (Central American Forearc) relative to southeastern Costa Rica (Panama microplate). In the backarc, this motion is expressed as extensional deformation, with strikes ranging from N140°E to N150°E, roughly parallel to the MAT azimuth in northwestern Costa Rica (N135°E). In contrast, the forearc accommodates deformation through lower-magnitude strains and E-W extension.