Deformation associated with sliver transport in Costa Rica: Seismic and geodetic observations of the July 2016 earthquake sequence.

Tectonic slivers form in the overriding plate in regions of oblique subduction. The inner boundaries of the sliver are often poorly deﬁned and can consist of well-deﬁned faults, rotating blocks or diffuse fault systems, which pass through or near the volcanic arc. The Guanacaste Volcanic Arc Sliver (GVAS) as deﬁned by Montero et al. , is a segment of the Central American Forearc Sliver, whose inner boundary is the ∼ 87-km-long Haciendas-Chiripa fault system (HCFS), which is located ∼ 10 km behind the volcanic arc and consists of strike slip faults and pull apart steps. We characterize the current ground motion on this boundary by combining earthquake locations and focal mechanisms of the 2016 Bijagua earthquake sequence, with the surface ground deformation obtained from Interferometric Synthetic Aperture Radar (InSAR) images from the ALOS-2 satellite. The coseismic stack of interferograms show ∼ 6 cm of displacement towards the line of sight of the satellite between the Ca˜no Negro fault and the Upala fault, indicating uplift or SE horizontal surface displacement. The largest recorded earthquake of the sequence was M w 5.4, and the observed deformation is one of the smallest earthquakes yet detected by InSAR in the Central American region. Forward and inverse models show the surface deformation can be partially explained by slip on a single fault, but it can be better explained by slip along two faults linked at depth. The best-ﬁtting model consists of 0.33 m of right lateral slip on the Ca˜no Negro fault and 0.35 m of reverse slip on the Upala fault, forming a positive ﬂower structure. As no reverse seismicity was recorded, we infer the slip on the Upala fault occurred aseismically. Observations of the Bijagua earthquake sequence suggests the forearc sliver boundary is a complex and diffuse fault system. There are localized zones of transpression and transtension and areas where there is no surface expression suggesting the fault system is not yet mature. Although aseismic slip is common on subduction interfaces and mature strike-slip faults, this is the ﬁrst study to document aseismic slip on a continental tectonic sliver boundary fault.


I N T RO D U C T I O N
When tectonic plates converge obliquely, the slip is partitioned into a trench perpendicular component taken up by subduction and a trench parallel component accommodated within the overriding plate (e.g. Jarrard 1986;McCaffrey 1996a). This leads to the formation of a tectonic sliver or microplate, which typically acts as a rigid block between the trench and a bounding fault system. In some cases, the boundary fault system consists of one or more strike slip faults, such as the Sumatra Fault that bounds the Sumatra forearc (e.g. Bradley et al. 2017), the Tarera-Aiduna fault system that bounds the Irian Jaya forearc (e.g. McCaffrey 1996b), and the Liquine-Ofqui faults that bounds the Central and Southern Chilean forearcs (e.g. Lavenu & Cembrano 1999). In other cases, the trench parallel component of motion is accommodated by block rotation, such as in the Mentawai forearc sliver off Sumatra (Berglar et al. 2017), the Aleutian forearc (Geist et al. 1988), the Cascadia forearc (e.g. Wells et al. 1998) and the segmented forearcs of El Salvador and Nicaragua (e.g. Turner et al. 2007;Alvarado et al. 2011). The boundary fault system often lies through the thermally-weakened volcanic arc, but in the case of the Guanacaste Volcanic Arc Sliver (GVAS) in Costa Rica, it is ∼10 km into the backarc (Montero et al. 2017).
Analysis of the earthquakes in forearc regions can be used to delineate fault systems and characterize the relative motion. Our analysis of the USGS earthquake catalogue shows that tectonic sliver boundaries commonly have moderate strike-slip earthquakes (up to M w 6.5) on mature boundary faults and smaller (<M w 5.0) earthquakes on the leading edges, trailing edges and on diffuse boundary fault systems (U.S. Geological Survey 2017). The pattern of seismicity typically shows strike slip events along the volcanic arc, with compression at the leading edge and extension at the trailing edge of the sliver (e.g. Wells & Coppersmith 1994;McCaffrey 1996a;Wang 1996;Lewis et al. 2008;Becker et al. 2010;Haq 2010). This is illustrated well by the Cascadia forearc, which experienced two reverse faulting events in Washington (the 2017 M w 4 Belfair earthquake and the 2014 M w 4.1 Seabeack earthquake) and a normal faulting earthquake in Oregon (the 1993 M w 6.0 earthquake, e.g. Braunmiller et al. 1995;Wang 1996). Aseismic motion is often observed on major strike slip faults [e.g. San Andreas Fault, Jolivet et al. (2015); Hayward Fault, Harris (2017); North Anatolian Fault, Rousset et al. (2016)]. However, no aseismic motion has yet been detected at tectonic sliver boundaries where slip rates are typically lower and aseismic motion may be harder to detect.
In this study, we focus on the GVAS Costa Rica, which formed as a result of the oblique plate convergence of the Cocos Plate beneath the Caribbean plate, at a rate of ∼88 mm yr −1 and angle of N23 • E (e.g. DeMets 2011, Fig. 1). The obliquity (angle between the strike of the plate boundary and the plate displacement vector) in this part of Costa Rica is <15 • (e.g. McCaffrey 1996a) and the GVAS moves as a relatively rigid block with a velocity of ∼8 to 11 mm yr −1 to the N45 • W, relative to the Caribbean Plate (e.g. Norabuena et al. 2004;DeMets 2011;Feng et al. 2012). The NE boundary of the GVAS is thought to be the Haciendas-Chiripa fault system (HCFS), which is a NW-SE striking fault system (Montero et al. 2017). The link between this fault system and the Central Costa Rica Deformed Belt (CCRDB, Marshall et al. 2000), is not well defined, due to the lack of surface expression and low seismicity rates. In this study, we used an interdisciplinary approach (seismology and geodesy) to investigate a seismic sequence that occurred on 3 July 2016 around the HCFS. Our goal is to analyse the deformation associated with the earthquake sequence using InSAR data, in order to obtain the fault source geometry and to relate the seismicity in the area to the mechanisms of sliver transport and regional tectonics. The largest earthquake in the sequence was M w 5.4, making this one of the smallest earthquakes ever detected with InSAR in the Central American region (e.g. Funning & Garcia 2018).

The HCFS
Based on GPS velocity vectors, Feng et al. (2012) proposed that the NE boundary of the GVAS is located within the volcanic arc. However, given the GPS station spacing in this area, the deformation could occur anywhere within a region spanning ∼30 km. More recently, Montero et al. (2017) suggested that the sliver motion is accommodated by the HCFS, which is located parallel to and ∼10 km NE of the volcanic arc. This backarc area consists of highly sheared Quaternary volcanic deposits and the trace of the HCFS offsets drainages and deflects streams by hundreds of meters, indicating recent lateral motion (Montero et al. 2017), (Table 1). The HCFS is composed of six NW striking faults, named from north to south: Haciendas, Caño Negro (CNF), Upala (UF), Chiquero, Cote-Arenal and Chiripa faults (Table 1 and Fig. 2a). The fault system includes two major steps, with pull apart basins between the CNF and Chiquero Faults (<100 m 2 ), and another between the Haciendas fault and the CNF. The 2016 Bijagua earthquake sequence occurred around the SE end of the CNF, near the step located between the Caño Negro and the Chiquero faults (see Fig. 2b). The Caño Negro and Chiquero faults are both right-lateral and mapped as near vertical with dips between 80 • and 90 • (Montero et al. 2017). The Upala fault is mapped as a South dipping reverse fault with an ignimbrite unit dipping >70 • S and with open cracks parallel to the fault trace.

Earthquake catalogue
The National Seismological Network (RSN) of Costa Rica has over 90 recording stations in the Guanacaste region and has detected persistent seismicity since the catalogue began in 1973 (Fig. 2a). 475 events were recorded and manually located from a minimum of eight stations and with a maximum distance to the station of ∼50 km from the hypocentre. The events were shallower than 20 km deep, with a minimum M w 1.6 and maximum of M w 5.4. For the location a Vp/Vs velocity ratio of 1.72 was used and the magnitudes were calculated from the instrument corrected displacement spectrum (using the flat spectral level and corner frequency) at six stations. The velocity model used is a 1-D P-wave velocity model of six layers generated for the Guanacaste Volcanic Arc, which is well constrained between ∼3 and 38 km depth (Araya et al. 2016). The top layer was difficult to constrain as it has inherent errors introduced by the assumption of laterally homogeneous crustal layers. The average relative errors of the relocated earthquakes with this velocity model are: ∼74 m east-west, ∼148 m north-south and ∼72 m in depth (Araya et al. 2016). From the catalogue of 475 manually located events in Araya et al. (2016), Montero et al. (2017) relocated 279 events between March 2000 and April 2018 within ∼6.0 km of the CNF using HypoDD (Waldhauser 2001).
In January 2002, a 4-d sequence of earthquakes occurred to the SW of the Caño Negro and Chiquero faults. The first and largest earthquake of the sequence was a M w 5.4 earthquake, (the 'main event'), followed by 150 detected aftershocks (Taylor et al. 2002). Taylor et al. (2002) reported that twenty-one of these aftershocks were felt by the local populations (e.g. in the towns of Upala and Bijagua), which had magnitudes between M w 3.0 and M w 4.4, and a depth range between 5.0 and 20 km (yellow events in Fig. 2b). Local residents reported surface cracking near the UF scarp (Montero et al. 2017). Unfortunately, no InSAR data is available to further investigate this sequence.
On 3 July 2016, the Bijagua earthquake sequence occurred SW of the HCFS, ∼7 km to the NE of Miravalles Volcano (Fig. 2b). On this day, at least 10 earthquakes had moment magnitudes between M w 2.5 and M w 5.4. The sequence occurred in an area ∼22 km long by ∼7 km wide, and between 1.8 and 8.4 km depth. This sequence is attributed to the HCFS because it aligns with the CNF, and has no apparent relation to volcanic activity. In particular, there is no seismicity recorded between this area and the Miravalles volcano summit and no temporal correlation between activity. The absolute errors for the individual manual location of the 2016 main event are 102 m N-S, 560 m E-W and 280 m in the vertical and after relocation using HypoDD (Waldhauser 2001) the relative errors for   those of the manual locations. Cross correlated waveforms between the source and a group of stations, were used to jointly calculate the distance, assuming the events travelled through the same path.
In this study, we focus on the two largest events of the 2016 earthquake sequence, for which focal mechanism solutions were obtained by Montero et al. (2017). The double couple earthquake focal mechanisms were determined by manually picking the polarities of the first arrival recorded at each station, using the location of the source and the incidence angles (Montero et al. 2017). We estimated additional source parameters using the empirical equations of Leonard (2010) for strike slip and dip slip faults. We used the moment magnitude (M w ), to calculate the seismic moment (M o ) (eq. 1) and to estimate the length (L), area (A) and slip (S) (eqs 2, 3 and 4, respectively, Table 2). All the other events in the sequence were of smaller magnitudes (<M w 4.0) and are unlikely to have contributed to the surface deformation. Event A (M w 5.4) and event B (M w 5.1) occurred ∼20 min apart in time, and <400 m apart horizontally, at depth of ∼1.9 and ∼8 km, respectively. Event (A) had a strike slip focal mechanism with near vertical nodal planes, striking NW-SE and NE-SW. The event (B) focal mechanism solution corresponds to a normal faulting earthquake with NS-striking nodal planes (Fig. 2b).

Processing and data
For this study, we processed SAR images of the Guanacaste region of Costa Rica acquired in 2016. The study area includes the HCFS, the 2016 Bijagua earthquake sequence, the Miravalles Volcano and national parks covered by tropical forests with dense vegetation. We used ALOS-2 PALSAR images from the Japanese Aerospace Exploration Agency (JAXA), which is equipped with an L band radar instrument (1.2 Hz sensor, ∼23.6 cm wavelength). The L-band wavelength penetrates vegetation and acquires more stable ground scatterers than shorter wavelength systems, hence it is less affected by decorrelation between images in densely vegetated regions. The major challenge with using this satellite is the long repeat intervals, and the fact that L-band is strongly affected by charged particles in the ionosphere. We also tested C-Band Sentinel 1A interferograms with a shorter time baseline (24 d), but they were not coherent enough to analyse (<30 per cent coherence), as is common between Latitude 15 • N and 15 • S, due to the dynamic atmosphere in these tropical latitudes and the dense vegetation and high water vapour (e.g. Funning & Garcia 2018).
We use six ALOS-2 images from the descending track 144, with a heading of ∼191 • and an angle of incidence ∼43 • . Repeat intervals were ∼112 to 196 d between acquisitions, at a sensing time of ∼17:26 UTM (i.e. Fig. 3a). The images used were taken between April 2015 and October 2016, from which we generated three interferograms before the seismic sequence, and three interferograms spanning the Bijagua earthquake sequence (Table 3), using the GAMMA remote sensing software package (Werner et al. 2000). To remove the topographic effect, we used a digital elevation model (DEM) with a 30 m resolution from the Shuttle Radar Topography Mission (SRTM, Farr et al. 2007, Fig. 3b). Given that we only have images from a descending track, the ground deformation measurements in this study are all in the descending line of sight (LOS) of the satellite (towards or away from the satellite), therefore we cannot obtain 3-D vectors for the ground displacements ).

Corrections
Each interferogram measures phase change between two consecutive images, which can be caused by a variety of factors: surface deformation, differences in the position of the satellite, thermal variations of the satellite during the acquisition, topographic effects, scatterer movement, and changes in atmospheric conditions (e.g. water vapour, aerosols; e.g. Bürgmann et al. 2000;Funning & Garcia 2018). To measure deformation associated with small earthquakes (<M w 6) from an interferogram, it is first necessary to reduce errors (e.g. Yu et al. 2018a,c). To reduce the white noise (thermal noise and loss of correlation), we filtered the interferograms with adaptive spectral filtering (Goldstein & Wernet 1998). Since coherence is particularly challenging in this densely vegetated region, we filtered twice with α = 0.4, before unwrapping the interferograms. We used a 30 per cent coherence threshold for masking pixels before unwrapping. We then unwrapped the phase using a minimum cost flow algorithm, with Delaunay triangulation (e.g. Costantini & Rosen 1999;Wegmuller et al. 2002, Fig. 3c). Atmospheric corrections can be particularly important in areas with high topographic gradients in tropical regions because low magnitude or deep seismic events can produce the same order of magnitude signal as tropospheric delays (e.g. Yu et al. 2018a). Water vapour concentration is vertically stratified through the atmosphere, where delays can be lower at summits and higher at low altitudes (e.g. Massonnet & Fiegl 1998). Thus, the pattern often correlates with the elevation contours. In contrast, if water vapour content is distributed in a turbulent way, it can create random, spatially correlated patterns (e.g. Hanssen 2001;Ebmeier et al. 2013;Parker et al. 2014).
To estimate the atmospheric contribution, it is necessary to characterize the pressure, temperature and water vapour at the acquisition time of each image. For this study, tropospheric delays were estimated using the Generic Atmospheric Correction Online Service for InSAR (GACOS, Yu et al. 2018b, see Fig. 3 D). GACOS uses the high-resolution weather model from the European Centre for Medium-Range Weather Forecasts (HRES-ECMWF) with a model spacing of ∼9 to 12 km every 6 hr and an Iterative Tropospheric Decomposition (ITD) model to separate the tropospheric turbulence and elevation dependent signals. It interpolates these using the SRTM DEM to produce maps of zenith delay at a resolution of 90 m (e.g. Yu et al. 2018b,c).
After the atmospheric corrections, we corrected for a longwavelength orbital error caused by the variations in the geometry of the satellite (position or small accelerations) between acquisitions and/or ionospheric effects (e.g. Fielding et al. 2018). We inverted for a two-dimensional linear empirical approximation for these variations (Parker et al. 2014) and subtract the results (Fig. 3e).
After the unwrapped interferogram is corrected for atmospheric and geometric effects, we masked out incoherent pixels using a threshold of <30 per cent coherence. To asses the overall level of coherence, we quantified the percentage of coherent pixels in our area of study, and found >90 per cent over the three coseismic interferograms (Table 3). In contrast, we obtained <27 per cent coherence in 48-d interferograms using Sentinel-1A images (C-band sensor, ∼5.6 cm wavelength). Table 2. Earthquake source parameters of the two biggest events from the 2016 Bijagua seismic sequence with area and slip estimated from empirical relations (see eqs 1-4) (Leonard 2010

Individual interferograms
We produced three interferograms (I1, I2, I3) spanning the 2016 earthquake sequence (Fig. 4 and Table 3). The wrapped interferograms are shown in Figs 4(a)-(c) and the unwrapped interferograms in Figs 4(d)-(f). A range decrease or displacement towards the satellite is shown in red, and blue represents a range increase or displacement away from the satellite.
The GACOS tropospheric delay maps used to correct the interferograms had maximum LOS displacement contributions of ∼2.8 cm in I-1, ∼8.7 cm in I-2 and ∼4.1 cm in I-3. The linear ramps had maximum contributions of ∼24, ∼7.8 and ∼3.4 cm, respectively over distances of <30 km. I-1 and I-2, have a clear linear ramp in a NE to SW direction, following the azimuth direction of the satellite. After corrections, the RMS for each interferogram decreased by 7.45, 2.51 and 0.22 cm, respectively, showing a greater decrease on the interferograms that showed larger RMS before correction (Table 3). There is no correlation between the coherence and the temporal baseline, suggesting seasonal variations are more significant than temporal baselines alone (e.g. Ebmeier et al. 2013).
All the obtained interferograms show ∼1 cycle or fringe (∼12 cm) of deformation in the line of sight (LOS) in an area of ∼30 km 2 between the CNF and the UF. The signal in I-2 is clearer than in I-1 or I-3; we attribute this to a smaller temporal baseline between acquisitions, even though I-1 and I-3 have better overall coherence. As I-2 has the shortest temporal baseline, with nearly the same magnitude and pattern signal as the other interferograms, we can infer that the bulk of the deformation occurred within the first ∼25 days after the Bijagua sequence, (i.e. from 3 July to 28 July). Interferograms I1 and I3, (Figs 4d and f), were generated from independent pairs of images, confirming the signal cannot be attributed to atmospheric artefacts. The peak deformation signal corresponds to a range decrease along the LOS, which implies a E-SE horizontal or uplift motion of the area between the Caño Negro and Upala faults, with the signal decreasing away from the faults.
Small differences in the deformation pattern can be seen between interferograms. For example, interferogram I3 has a small negative range change (red) on the SW side of the CNF, that is not present in I1 or I2. As I2 and I3 share the image acquired on the 7th of April 2016, we conclude that this signal is caused by a turbulent atmospheric effect on the 20th of October 2016 (Fig. 4f).

Stacked data
To further reduce atmospheric noise and other random noise, we averaged the interferograms assuming the same deformation signal is present in each (Table 3, e.g. Emardson et al. 2003). We generated two stacks: (1) a 1 yr stack of three interferograms before the 2016 Bijagua earthquake sequence (from April 2015 to April 2016) and (2) a stack composed of three interferograms that span the earthquake sequence (from April 2016 to October 2016) (Figs 5a and  b, respectively). The pre-seismic stack shows no signal of ground deformation in the area. The coseismic stack has a deformation pattern similar to the three individual interferograms. The coseismic stack better defines the area of deformation than the individual coseismic interferograms, and shows that the deformation signal is concentrated between the UF and CNF, but extends towards the Chiquero Fault (Fig. 5b). On the SW side of the CNF, the small turbulent signal observed in I3 remains, but with a smaller magnitude, reinforcing the conclusion that the signal is caused by turbulent atmospheric effect. On this side of the CNF we observe subsidence or NW horizontal motion, that extends along the CNF and towards the Chiquero fault (Fig. 5b).
The RMS of the coseismic stack is <1 cm, which is significantly smaller than the maximum deformation signal (∼6 cm, Table 3). In contrast, the RMS noise of the individual interferograms before the corrections were between 3 and 9 cm, meaning the deformation pattern would have been hard to distinguish.

S O U RC E M E C H A N I S M
We investigate the source mechanism for the 2016 Bijagua earthquake sequence by producing forward models of the ground deformation expected based on the seismic catalogue and inverse models using known fault geometry. Our aim is to compare the geodetic models to the seismic and geologic information to characterize and constrain the parameters of the active faults and their relationship to the HCFS and the regional tectonics. First, we generate forward models using the earthquake locations, focal mechanisms, and magnitudes to estimate the expected deformation for the two biggest events (A, B). Then we inverted the InSAR data using a linear inversion method to obtain the amount and distribution of slip on the known faults.

Forward modelling
We estimate 3-D ground displacements for events A and B (1-2) using the Okada (1985) model for shear motion on a rectangular plane dislocation in a homogeneous elastic half-space, using Lamés constants λ and μ of 3.2 GPa. Then, we projected the 3-D surface displacement field into the satellite LOS to produce a synthetic interferogram using a heading of −169 • and an incidence angle of 43 • . The nine source parameters used are latitude and longitude, length, strike, dip, rake, slip, top and bottom depth ( Table 2). The source parameters were obtained from seismic analysis of the earthquakes (focal mechanisms) and earthquake scaling relations (slip and length) assuming a square rupture area (eqs 1-4). Due to the high density of seismic stations in the area, the relative locations are more accurate than located with the global network (see Section 1.1). For the main event (A), we assumed a square rupture patch of 25.5 km 2 and a horizontal slip of 0.19 m obtained from scaling relations based on Leonard (2010) (see Table 2). This is only a small section (5 km length) of the much larger CNF. The forward model is based on the seismic location of the main event and errors in the absolute seismic location and possibly the dip of the fault plane mean that the deformation is not aligned with the surface trace.
The forward model of the largest earthquake (A), shows four main lobes of deformation near the SE end of the CNF trace, as is characteristic for a N-S strike slip fault (Fig. 6a). The deformation lobes on the SW block have diffuse patterns, that decrease gradually over a distance of ∼7.0 km. The lobes over the northern block have sharper deformation patterns close to the fault trace. The southern and eastern lobes show a range decrease of ∼5.4 cm and the northern and western lobes show a range increase of ∼7.3 cm. The difference in the signal pattern between the SW and NE blocks from the CNF is influenced by the sensitivity of the LOS to the vertical and horizontal components of displacement, creating an asymmetric pattern .
The focal mechanism for event B (M w 5.1) shows a normal fault striking NS, with one of the possible fault planes dipping to the west (B1) and a second dipping to the east (B2), consistent with the pull apart structure described in section 1.1. We made forward models for both nodal planes, (Figs 6B1 and B2). Fig. 6(B1), shows the forward model for the west dipping nodal plane solution with a maximum ground displacement of ∼1.4 cm away from the LOS. This ground displacement is located on the CNF, and covers an area of 13 km × 9 km (115 km 2 ). Fig. 6(B2), shows the forward model for the East dipping nodal plane, which shows a similar pattern of deformation as in B1 but the deformation is now located between the UF and Chiquero faults. The predicted displacements are much smaller (<1 cm) than for event A due to the much greater depth.
To illustrate the deformation pattern expected from the entire sequence, we combined the displacements of the individual events: C1 represents the combination of W-dipping B1 with event A, and C2 representing the displacements for the E-dipping B2 with event A (Figs 6C1 and C2). The pattern of LOS deformation predicted by the forward models (Fig. 6) is very different to the observed pattern in the coseismic stack (Fig. 5b). Although the predicted displacement is ∼6 cm of uplift, the spatial pattern is very different, as it covers a much smaller area than the deformation observed and is confined to just one block. Event A is shallower and dominates the predicted deformation pattern, with uplift or SE horizontal displacement at the SE end of CNF. One possibility is that aseismic slip (either coseismic, after slip or a slow slip event) contributed to the deformation pattern but was not recorded by the seismic network.

Inverse models
Next, we use a linear inversion procedure to estimate the amount and distribution of slip along faults of known geometry. We used the MATLAB coded 'SlipInv' (Funning 2005) that estimates the distribution of slip using a linear inversion of the geodetic data. We assumed the dislocation occurred along rectangular fault planes, on a homogeneous elastic medium, with the elastic Lamés constants λ, and μ, used for Section 4.1. We use the geometry of Montero et al. (2017) to define the fault geometry to test models that consider slip on one fault (CNF or UF), or on both faults. We discretize each fault plane into patches of 200 × 200 m (Fig. 7a, Table 4). The ground displacement data used for the inversion is the coseismic stack (Fig. 7a) subsampled into 1521 points, from the 66524 interferogram points, using a quadtree algorithm (Jonsson et al. 2002).
The fault geometry of the CNF is considered vertical (90 • dip), with a strike of S50 • E, right lateral strike slip motion, a length of 10 km and bottom depth of 4.8 km (Fig. 7b). There are no focal mechanisms or recorded seismicity near the Upala fault and so the dip at depth is unknown. We did a grid search for dip angle from 20 • to 70 • in steps of 5 • . For both models that consider the UF, the lowest RMS is obtained for a 35 • dip of the UF. The UF geometry used is 9 km length, striking ∼S68 • E. For this geometry, the CNF and UF intersect at a depth of ∼2.6 km (Fig. 7c). We test models for (a) pure right lateral slip on the CNF only, (b) pure thrust motion one the UF only and (c) slip on both the CNF and UF. We also did a grid search for the rake of both faults, from 160 • to 200 • for the CNF and from 70 • to 110 • for the UF and find that the slip on the CNF is near pure right lateral slip and near pure reverse slip on the UF, so the rake was fixed to exclude oblique motion in the subsequent models.
The inversion results show the key features of the deformation pattern can be explained by either single fault models with a similar misfit (∼1.64 cm), but that each model explains different aspects of the deformation pattern. The Caño Negro single fault model has two slip patches, one main slip patch has a rupture area of 6.6 km 2 at ∼3.0 km depth, and has maximum 0.71 m of right lateral slip located under the main deformation lobe between the CNF and the UF. The second patch has a rupture area of 0.66 km 2 at ∼0.80 km depth, located under the NW deformation lobe (Figs 7b, e and 8a). The Upala single fault model shows an elongated slip patch of 6.2 km 2 at ∼ 2.0 km depth and a maximum of 0.45 m of reverse slip located at the bottom of the fault plane, near where the Upala fault meets the CNF (Figs 7c, f and 8b). This model explains the deformation pattern between the faults better than the Caño Negro single fault model, but it does not explain the second lobe of deformation located to the W of the UF (Fig. 7f).
The two-fault model explains both deformation lobes and with a smaller misfit (1.47 cm) than the two single fault models (Table 4, Figs 7d, g and 8c). The best fit result has slip on the Caño Negro Fault distributed in two patches on the fault plane, with a right lateral slip direction (Fig. 8). The SE patch is 5.9 km 2 with a maximum slip of ∼0.19 m at ∼4 km depth, equivalent to 3.65 × 10 16 Nm geodetic moment or M w 4.9. The NW patch is ∼4 km 2 , with a maximum ∼0.33 m slip at ∼1 km depth, equivalent to 4.22 × 10 16 Nm geodetic moment or M w 4.9. The total geodetic moment on the CNF is 7.88 × 10 16 Nm, equivalent to a M w 5.1 earthquake. For the Upala fault, the two-fault model shows one elongated patch of reverse slip, this type of motion has not been observed at any recorded seismic event near this area. The slip patch area is ∼6.2 km 2 with a maximum of ∼0.35 m of reverse slip at a depth of ∼2.2 km. The total geodetic moment associated with the Upala fault is 6.9 × 10 16 equivalent to a M w 5.0 earthquake.

D I S C U S S I O N
The seismic records are used to calculate the ground deformation expected from the recorded earthquakes. In the study region, we know there are higher uncertainties from events located at depths <3 km, when using the 1-D velocity model from Araya et al. (2016). As the main event is a M w 5.4 earthquake, located at 1.9 ± 0.28 km depth, ground deformation is expected. The deformation observed with geodesy can be used to locate the source with lower uncertainties and better characterize the source. The second largest event is a M w 5.1 earthquake located at 8 ± 0.55 km depth, at a depth better constrained by the 1-D velocity model, from which very little (<1.5 cm) ground deformation is expected. The coseismic interferograms spanning the 2016 Bijagua earthquake sequence show ∼6 cm between the Caño Negro, Upala and Chiquero faults (Fig. 9b), in the line of sight of the descending track of the satellite, which could be uplift or SE horizontal motion. Here we consider the relationship between the deformation and the two largest earthquakes of the seismic sequence (events A and B) that occurred at the SE of the CNF trace, as they are the only ones large and shallow enough to generate detectable ground deformation (Funning & Garcia 2018). In Fig. 9, we plot the displacement in the LOS as a 13-km-long cross section (P' -P") from the NW-SE striking Caño Negro towards the Upala Fault, showing the hypocentre of the seismic events. The cross section shows clearly that the location of deformation coincides with the area between the CNF and UF. From the geodetic inversion, we obtained a maximum reverse slip of 0.35 m on the UF, which is slightly larger than the slip on the CNF (0.33 m). The profile also shows that significant change of the displacement coincides with the fault traces, with the peak displacement between the two faults. These means there probably was no surface rupture and that the source is deep, Fig. 9(b).
The distribution of the slip in the best-fitting geodetic model shows right lateral slip on the CNF, consistent with focal mechanism of Event A. The total geodetic moment between the two faults Downloaded from https://academic.oup.com/gji/article-abstract/220/1/585/5601383 by University Library user on 18 November 2019  is 1.48 × 10 17 Nm, equivalent to a M w 5.3 event, which is comparable but slightly smaller than the seismic estimate for event A (1.58 × 10 17 Nm). Event B was deep (∼8 km) and probably did not cause any detectable deformation. The direction of slip modelled on the Upala fault is reverse and is not consistent with any recorded earthquakes during the 2016 Bijagua sequence, for this reason we assume the slip occurred aseismically. Although, our satellite observations cannot pinpoint the precise timing of the aseismic slip, there is no recognizable deformation signal on the pre-seismic stack, and as the observed deformation appears on the shortest time baseline interferogram (I2, Figs 4b and e), we assume it was triggered by the strike slip event on the CNF and occurred sometime in next the 25 d (Fig. 9a, red star).
The seismic location of Event A lies ∼2.6 km from our geodetic solution towards N31 • W, and ∼300 m deeper. Although the seismic relocation process reduces the relative errors between events in the sequence, there is often still a shift in absolute locations. However, the difference in location could arise because the seismic location is of the first motion of the rupture, and the geodetic location is the centroid of the rupture, in which case, the difference between locations indicates the direction of the rupture during the event. For small events like these, it is more likely that the difference between the locations are due to errors in the seismic locations because 1-D velocity models do not consider lateral heterogeneities, and in particular for the model used, the upper 3 km are not well constrained (Araya et al. 2016). As the in situ location given by InSAR is more accurate (Weston et al. 2011), we relocate the entire  (Fig. 9). sequence based on the more accurate geodetic location for Event A (Fig. 9) and using the relative location of the smaller events from the Hypo DD relocation. In Fig. 9, we show the new geodetic locations (dark read circles) of the Bijagua tectonic sequence based on the difference (red line) between the double difference location and the geodetic location of the main event.
The distance on the surface between the CNF and UF is ∼3 km, and based on our slip inversions, the slip on the Upala fault occurred at the depth near the intersection with the Caño Negro fault. The relationship between these faults can be explained by a positive flower structure consistent with a transpressional stress field. Aseismic displacement on a linked thrust could be triggered by motion on a strike slip fault, if there is low friction on the reverse fault plane, or if the reverse fault was already close to failure (e.g. Bayasgalan et al. 1999;Fielding et al. 2004). This fault geometry is similar to that of the 1998 M w 6.6 Fandoqa strike-slip earthquake, Iran, which also had associated post-seismic afterslip on a linked thrust (e.g. Fielding et al. 2004). In both cases, a strike slip earthquake triggered aseismic motion on the base of an intersecting thrust fault, within a positive flower structure. Our observations suggest the Bijagua seismic sequence occurred in a local transpressional zone, but as a M w 5.1 normal motion earthquake (event B) also occurred during the seismic sequence, it is likely that there are also local transtensional zones, consistent with the observed pull apart basin.

C O N C L U S I O N S
We studied the coseismic deformation caused by a 2016 seismic sequence on the HCFS associated with the GVAS boundary. This analysis was performed using ALOS-2 satellite InSAR images over a period of 280 days, to compare the results with the seismologic and geological data. Although there are challenges to using InSAR in Costa Rica (at a tropical latitude and with dense vegetation), we have detected deformation patterns on interferograms of an earthquake sequence with a M w 5.4 main event, demonstrating the potential of the technique in the area. This is one of the smallest earthquakes ever to have been detected with InSAR in the Central American region. However, the low repeat frequency and availability of ALOS-2 still presents a challenge in the region. In contrast, Sentinel 1A has a better repeat frequency but suffers from decorrelation between Lat 15 • N to 15 • S (e.g. Funning & Garcia 2018).
In the past 20 yr, seismicity along the HCFS has been concentrated on the SW side of the Caño Negro and Chiquero faults, mainly with strike slip focal mechanisms, but also with moderate normal motion focal mechanisms. The main deformation observed on three independent interferograms from 2016 is between the Caño Negro fault (CNF) and the Upala Fault (UF), with a maximum ∼6 cm of uplift or SE horizontal displacement. The forward model deformation pattern based on seismic scaling relations and the focal Figure 9. (a) Summary map of the Bijagua seismic sequence. The red star represents the seismic location and the blue star is the geodetic centroid of the deformation for the main event. In dark red points are the new locations for the rest of the seismic sequence based on the new geodetic location for the main event. The focal mechanism (a) is the main earthquake M w 5.4 from 2016 earthquake sequence and (b) the M w 5.1 earthquake from the same sequence. (b) Cross section between P'-P", in blue and black lines show the displacement in the line of sight of the satellite of the coseismic stack and the topography, respectively. The red line shows the model of displacement in the LOS of the two-fault model. The grey band shown in the profile is the RMS range of the phase change from the coseismic stack. (c) Shows the inferred motions of the CNF and UF on the bottom of the cross section and the new locations for the earthquakes of the 2016 Bijagua earthquake sequence based on the geodetic analysis of the main event, and the focal mechanisms for events A and B. mechanisms, did not fit the observed deformation pattern. The results from an inversion of the slip distribution on the known fault geometry found a maximum right lateral slip of 0.33 m on the CNF plane and 0.35 m of reverse slip in the UF. We were not able to associate the motion on the UF with any recorded seismic event, even though the geodetic moment on this fault represents 37 per cent of the total. For this reason, we concluded the slip on the UF occurred aseismically but triggered by the earthquake sequence. This demonstrates the benefit of combining geodetic and seismic data. This study is the first to report aseismic slip on a sliver boundary, this is relevant because aseismic slip is typically associated with major transform faults or subduction interfaces.
We propose that the Caño Negro and Upala faults are linked through a developing positive flower structure associated with the sliver boundary, which accommodates the trench parallel slip from the oblique subduction. This flower structure may repeat or continue parallel to the Caño Negro fault, but as no seismicity or deformation has yet occurred to prove its existence, we limit our mapping to a length of ∼8 km for now. The flower structure inferred from this study has a geometry similar to that found in the fold-and-thrust belt of Iran, which also displays aseismic slip on splay reverse faults linked to a major strike slip fault (Fielding et al. 2004).
Through this study we observed the GVAS has a complex sliver boundary, which has primary NW right lateral striking faults with a pull apart basin between the steps of the faults and parallel reverse splay faults linked at depth. These observations demonstrate the complexity of this inner sliver boundary, which has active transtension and transpression occurring on different local sections of the same fault system. We recommend further geophysical studies to improve the knowledge of the structure of the HCFS and the associated seismic hazards.

A C K N O W L E D G E M E N T S 6 A c k n o w l e a d m e n t s
We would like to thank the University of Costa Rica and the National Seismological Network, for allowing the use of their earthquake catalogue. We want to thank CEOS Disaster Risk Reduction Pilot Project for the ALOS-2 satellite images from the Japanese Aerospace Exploration Agency. We thank Gareth Funning for Slip-Inv MATLAB code. We want to thank Luke Wedmore for his suggestion for the earthquake scaling relations equations. We will also like to acknowledge the University of Costa Rica (UCR) and the NERC Centre for the observation and Modelling of Earthquakes, Volcanoes and Tectonics (COMET) for the financial support. For the interferogram processing we used GAMMA (Werner et al. 2000) and for the stacking and figures we used MATLAB (MathWorks 2005) and GMT (Wessel & Smith 1998).