Geophysical Journal International

SummarySeismic anisotropy is key to constrain mantle flow, but it is challenging to image and interpret it. Existing large-scale tomography models of seismic anisotropy typically show large discrepancies, which can lead to completely distinct geodynamical interpretations. To better quantify the robustness of anisotropy tomography, we create a 2-D ridge-to-slab geodynamic model and compute the associated fabrics. Using the resulting 21 elastic constants we compute seismic full waveforms, which are inverted for isotropic and radially anisotropic structure. We test the effects of different data coverage and levels of regularisation on the resulting images and on their geodynamical interpretation. Within the context of our specific imposed conditions and source-receiver configuration, the retrieved isotropic images exhibit substantial artificial slab thickening and loss of the slab’s high velocity signature below ∼100 km depth. Our results also show that the first order features of radial anisotropy are well retrieved despite strong azimuthal anisotropy (up to 2.7 %) in the input model. On the other hand, regularisation and data coverage strongly control the detailed characteristics of the retrieved anisotropy, notably the depth-age dependency of anisotropy, leading to an artificial flat depth-age trend shown in existing anisotropy tomography models. Greater data coverage and additional complementary data types are needed to improve the resolution of (an)isotropic tomography models.

SummaryIn the context of Bayesian inversion, we consider sequential Monte Carlo (SMC) methods that provide an approximation of the posterior probability density function and the evidence (marginal likelihood). These particle approaches build a sequence of importance sampling steps between gradually-tempered distributions evolving from the prior to the posterior PDF. To automate the definition of the tempering schedule, adaptive sequential Monte Carlo (ASMC) allows tuning the temperature increments on-the-go. One general challenge in Bayesian inversions is the computational burden associated with expensive, high-fidelity forward solvers. Lower-fidelity surrogate models are interesting in this context as they can emulate the response of expensive forward solvers at a fraction of their cost. We consider surrogate modeling within ASMC and introduce first an approach involving surrogate modeling only, in which either prior samples are used to train the surrogate, or the surrogate model is re-trained by updating the training set during the inversion. In our implementation, we rely on polynomial chaos expansions for surrogate modeling, principal component analysis for model parametrization and a ground-penetrating radar cross-hole tomography problem with either an eikonal or finite-difference time-domain solver as high-fidelity solver. We find that the method based on re-training the surrogate during the inversion outperforms the results obtained when only considering prior samples. We then introduce a computationally more expensive multifidelity approach including a transition to the high-fidelity forward solver at the end of the surrogate-based ASMC run leading to even more accurate results. Both methods result in speed-ups that are larger than one order of magnitude compared to standard high-fidelity ASMC inversion.

SummaryDense array observation and seismic interferometry have revolutionized the imaging schemes of the earth structure. It is becoming possible to directly obtain the lateral variation of the earth's structure by applying array-based methods such as the cross-correlation beamforming (CBF) of the ambient noise to the subsets of the dense array, without tomography. CBF has been proven to extract the azimuth-averaged multi-mode surface wave dispersion curves. However, the resolution of the dispersion image generated by conventional CBF is low at high frequencies in the frequency-velocity (f-v) domain. Moreover, the irregular array geometry and uneven source distribution would bias the result of CBF, especially for the estimation of azimuth-dependent velocity. In this paper, two beamforming schemes are suggested to improve the resolution of multi-mode dispersion images in the f-v domain. Firstly, the geometrical spreading of the wavefield is corrected to enhance the amplitude at high frequency (or large distance) and thereby improve the resolution of the dispersion image at high frequency. We call this scheme weighted correlation beamforming (WCBF). The azimuth-averaged velocity can be estimated with sufficient resolution using WCBF by stacking the beamforming output at each azimuth. We show that WCBF is the 2D Fourier transform of the spatial wavefield from the viewpoint of the wavefield transform. Secondly, a modified beamforming scheme (MCBF) is suggested to reduce the effect of uneven source and/or irregular array geometry. The delay and summation in MCBF are performed only for plane waves incident from the stationary phase region. The azimuth-dependent velocity can therefore be estimated by MCBF with less dependence on the array geometry, as well as on the uneven source distribution. In terms of the estimation of azimuth-averaged phase velocity, we show the F-J method, another array-based method for extracting multi-mode surface waves from ambient noise using the Fourier-Bessel transform, is the azimuth-averaged version of WCBF. The reliability of WCBF and MCBF is verified based on the synthetic and field data using the array with different geometry. The dispersion image of multi-mode Rayleigh wave phase velocity at local and regional scales can be generated by WCBF or MCBF with high resolution. In particular, multi-mode dispersion curves at the local scale can be measured by MCBF with sufficient accuracy using quite short recordings from hours to days. This offers the possibility of a rapid assessment of the media properties.

SummaryTo understand the seismic hazard of a subduction zone, it is necessary to know the geometry, location, and mechanical characteristics of the interplate boundary below which an oceanic plate is thrust downward. By considering the azimuthal dependence of converted P-to-S (Ps) amplitudes in receiver functions (RFs), we have detected the interplate boundary in the Makran subduction zone, revealing significant seismic anisotropy at the base of the accretionary wedge above the slab before it bends down beneath the Jaz Murian basin. This anisotropic feature aligns with a zone of reduced seismic velocity and a high primary/secondary wave velocity ratio (Vp/Vs), as documented in previous studies. The presence of this low-velocity highly anisotropic layer at the base of the accretionary wedge, likely representing a low-strength shear zone, could possibly explain the unusually wide accretionary wedge in Makran. Additionally, it may impact the location and width of the locked zone along the interplate boundary.

SummaryA major earthquake shook the Chinese county of Maduo, located in the Songpan-Ganzi terrane on the Tibetan Plateau, on 21 May 2021. Here, we investigate the postseismic deformation process of this event, with the aim to understand the fault geometry, friction behaviour and regional rheology. To keep the self-consistency between co- and post-seismic deformation models, we first constrain the fault geometry and coseismic slip model of this event, which are directly used in modelling the postseismic deformation. The coseimsic slip model reveals that the majority of coseismic slip is confined at the middle (3-15 km) of the brittle layer, leading to significant shallow slip deficit. Secondly, we obtain the postseismic deformation in the first 450 days following the 2021 Maduo earthquake using the GPS and InSAR displacement time series data. Thirdly, a combined model incorporating afterslip and viscoelastic relaxation is built to explain the observed postseismic deformation. Our results suggest that the viscoelastic relaxation effect should be considered in the observation period, in order to avoid the unphysical deep afterslip in the ductile lower crustal layer. Combined analysis on viscosities inferred from this study and previous studies suggests a weak lower crust with steady-state viscosity of 1018 ∼ 1019 Pa s beneath the Songpan-Ganzi terrane, which may give rise to the distributed shear deformation and the development of subparallel secondary faults within the terrane. Besides, the inferred afterslip on uppermost patches of the middle fault segment suggests a rate-strengthening frictional behaviour that may be related to the coseismic slip deficit and rupture arrest of the Maduo earthquake.

SummaryAmbient noise technology can efficiently extract surface wave signals from seismic background noise and has been extensively utilized in imaging lithospheric structures. However, retrieving crustal body wave signals, such as PmP or SmS phases, still poses a challenge. Only a limited number of reports have successfully extracted these regional-scale body wave signals from ambient noise in only a few limited study areas. It remains unclear why these signals are difficult to retrieve from ambient noise data.To investigate the mechanism of recovering body wave signals in noise cross-correlations, we calculate cross-correlation functions at four regions and observe the similarity of the recovered body waves. Through a series of synthetic simulations, we demonstrate that the appearance of body wave signals in noise cross-correlations is closely related to the distribution of noise sources. Among these signals, the post-critical SmS wave proves to be the most readily recoverable from ambient noise data, primarily stemming from distant sources. In contrast, the recovery of P-wave requires the array to be in proximity to the sources. Our experiments also reveal that the main origin of PL waves is the multiple reflections of S-waves propagating in the crust.

SummaryFor a more quantitative discussion of slow earthquake activity, we evaluated the detectable limits of very low frequency earthquakes (VLFEs), which are seismic slow earthquakes observed in very low frequency (< 0.05 Hz) bands in the Nankai subduction zone. We performed numerical simulations using a local three-dimensional model and used the observed noise level of permanent broadband seismometers. First, we investigated the effects of the source-time functions on the maximum amplitudes of the VLFE signals at a certain station. The maximum amplitudes of the VLFE signals were controlled by the VLFE moment rate. The detectable limit of VLFEs at each source location can be defined as the lowest moment rate of detectable VLFEs, which radiate signals larger than the noise levels of any component at ≥ 3 stations. For inland seismometers only, the detectable limits of VLFEs at deep (30–40 km) and shallow (≤ 10 km) depths were 1012–1012.3 and 1012.7 Nm/s, respectively. Due to the geometrical spreading of VLFE signals and large noise levels in horizontal components, offshore seismometers improved the detectability of shallow VLFEs in regions where seismometers were densely deployed. Based on our detectability and published catalogs, shallow slow earthquakes are less active south-southwest off the Kii Peninsula, where geodetic studies expect mechanical coupling.

SummaryWe have derived an analytical approximate expression to estimate the delay in diving seismic waves due to thin layers of CO2. The expression is valid for high frequencies and can be used to estimate the delay in diving waves at seismic frequencies for large separations between the source and receiver (offset). The approximation may be used to assess CO2 detection limits using diving waves and to support survey planning for CO2 monitoring and full-waveform inversion (FWI) cycle skipping analysis. In this study, we analyze the diving-wave response to a thin layer of CO2 for band-limited data using acoustic finite-difference modeling, and compare the results against the analytical calculations. We find that the responses are offset-dependent and related to double- and single-leg interactions between the diving waves and the CO2. To test the methods, we created a synthetic representation of the 2010 subsurface conditions for the top CO2 layer at the Sleipner storage complex in the North Sea, by combining base and monitor post-stack seismic data with field velocity trends. Using the acoustic finite-difference method, we model pre-stack data that captures the complexity of field data and demonstrate the use of the diving-wave delay for CO2 migration monitoring and CO2 thin layer detection.

SummaryMagnetotelluric data from Mount Tongariro has been analysed using an unstructured tetrahedral finite-element inversion code that incorporates topography, which was not included in previous analysis of these data. Incorporating topography adds information, which stabilises the resistivity inversion modelling, and for the first time allows details of the shallow hydrothermal system and its relationship with the underlying magmatic system to be resolved. Specifically, an electrically conductive zone between 4 and 12.5 km depth marks the underlying magmatic system, which is shown to directly connect via conductive pathways to the area where the most recent phreatic eruptions at Tongariro occurred in 2012. The resultant phreatic eruptions in August and November 2012 showed no new magmatic component to the eruption deposits. Nevertheless, by combining the magnetotelluric resistivity image with relocated seismicity, we can see that seismicity (a proxy for magma ascent) migrated from the top of the magmatic system into the hydrothermal system in the months preceding these eruptions. Magmatic interaction with the extant hydrothermal system likely caused the over pressurisation for the phreatic eruption. This work highlights the utility of combining geophysical methods and petrological data to constrain phreatic eruption processes.

SummaryHydraulic Fracturing (HF) often stimulates the local earthquake productivity which provides a unique opportunity to characterize the crustal heterogeneities, reservoir properties, and fluid injection effects. However, the velocity models acquired solely based on the arrival time records are often undermined due to the seismic network coverage and interpolation techniques. Instead, we adopt the waveform-based approach to apprehend; (1) structural heterogeneities, (2) reservoir distribution, and (3) signatures of the injected fluid in the Weiyuan shale gas field (WSGF). We categorize the waveforms into dominant high and low frequencies based on the qualitative inspection and frequency index analysis of the seismic waveforms. We first inspect the waveform to access the potential controlling mechanisms (source, site, and path effects) at both single and multiple stations in different azimuthal orientations. As a result, we find the path effect as a dominant factor to influence the waveform characteristics, e.g., S-wave amplitude, and frequency. Subsequently, to localize the path effect, we conduct an in-depth examination of events within 10 km of each seismic station and classify the waveform records using their frequency indices. Notably, certain stations record a significant proportion of low-frequency waveforms (LFWs) (up to 20%), while others have limited occurrences (∼1%) indicating suspected anomalous zones. Afterward, we identify two suspected anomalous zones based on LFWs intensity and ray tracing map. Both zones are in close proximity to fault zones and preserved reservoirs with no HF activities, where fault damage zones or the fluid-rich reservoir may contribute to our observed LFWs.

SummaryThe south-eastern offshore of the Boso peninsula in Japan periodically experiences short-term slow slip events (SSEs) every few years. On June 2018, a SSE occurred with the maximum surface horizontal displacement reaching up to 4.7 cm by according to the operational global navigation satellite system (GNSS) network. This study performed a time series analysis of interferometric synthetic aperture radar (InSAR) with Sentinel-1 SAR images to investigate detailed spatial pattern of surface displacements caused by the SSE. With the assistance of an atmospheric delay correction with a regional numerical weather model output, the InSAR time series analysis successfully captured displacement signals in three paths, whose maximum amplitudes in line-of-sight directions were 1.46 cm, 1.86 cm and -0.80 cm. A checkerboard test revealed that the resolution of the slip inversion was higher when InSAR was used than that using GNSS, especially in and around the inland. The slip inversion with the actual displacement data derived from the InSAR time series analysis was performed with the L-curve optimization, showing that the estimated slip area was concentrated offshore south-eastward from the Boso peninsula with the maximum slip of 5 cm and the estimated moment magnitude of 6.4. As similar to previous SSEs in the Boso peninsula, a seismic swarm simultaneously occurred in the down-dip area adjacent to the estimated slip with the SSE occurrence, suggesting a different friction characteristics between them. This study demonstrates usefulnesses of the InSAR observation for capturing detailed spatial characteristics of small-displacement events like SSEs and of the hybrid use of the externally derived delay correction with the time series analysis to improve the displacement detection accuracy.

SummaryLow-frequency earthquakes (LFEs) are seismic phenomena with the shortest timescale among various slow earthquakes observed on broadband time scales. To understand the nature of such a broadband slow phenomenon, it is important to investigate the rupture evolution process of individual slow events, such as LFEs. Here, we investigated the moment-duration relationship of LFEs at plate interfaces and volcanic regions, and showed that the moment-duration relationship of both tectonic and volcanic LFEs is characterised by a moment proportional to the cubic duration, similar to that in ordinary earthquakes. The difference between our obtained moment-duration relationship and the broadband scaling suggests that the evolution process of LFEs may not be controlled, but only triggered by the slow earthquakes with longer durations, such as slow slip events driven by aseismic diffusion. The seismic moments of the LFEs are approximately three orders of magnitude smaller than those of ordinary earthquakes with similar durations. This result indicates that LFEs have rupture growth similar to that of ordinary earthquakes, although the rupture velocity and/or stress drop are much smaller. Considering the hypocentre spread of LFEs, the estimated rupture velocity and stress drop were approximately 100 m/s–1 km/s and 2 kPa–1 MPa, respectively. Additionally, the estimated moment magnitudes are much larger than the local magnitudes determined based on the maximum amplitudes, which is due to the longer durations and resultant smaller amplitudes of LFEs than those of ordinary earthquakes.

SUMMARYWe present a time-domain distributional finite-difference scheme based on the Lebedev staggered grid for the numerical simulation of wave propagation in acoustic and elastic media. The central aspect of the proposed method is the representation of the stresses and displacements with different sets of B-splines functions organized according to the staggered grid. The distributional finite-difference approach allows domain-decomposition, heterogeneity of the medium, curvilinear mesh, anisotropy, non-conformal interfaces, discontinuous grid, and fluid-solid interfaces. Numerical examples show that the proposed scheme is suitable to model wave propagation through the Earth, where sharp interfaces separate large, relatively homogeneous layers. A few domains or elements are sufficient to represent the Earth’s internal structure without relying on advanced meshing techniques. We compare seismograms obtained with the proposed scheme and the spectral element method, and we show that our approach offers superior accuracy, reduced memory usage, and comparable efficiency.

SUMMARYCorrecting atmospheric effects in measurements of vertical acceleration is challenging, particularly at frequencies above 0.3 mHz. Corrections based on four-dimensional atmosphere models should be the most effective, but are limited in temporal and spatial resolution. So these models are commonly combined with a constant admittance between air pressure and gravity, though the correction can be further improved by allowing a frequency-dependent admittance. This paper studies the effectiveness of such procedures given that air pressure admittance varies with time because of variability in the underlying local atmospheric process, and finds that admittances estimated by cross-spectral analysis gives better results than a constant admittance does. The best results are achieved by combining a 4D model with admittances estimated from time series that cover the event of interest.

SUMMARYDike propagation is a mechanism for more rapid ascent of felsic magmas through the crust than is possible via diapirs or percolative flow. As it ascends, the magma undergoes complex physical and chemical transformations induced by decompression and cooling. These processes dramatically change the magma density and viscosity, which in turn affect magma ascent rate and the depth at which the dike arrests. We present a mathematical model of dike propagation for silicic magmas taking into account the presence of multiple volatile species (H2O and CO2), bubble growth, heat advection and loss, crystallization and latent heat release. We consider conditions for dikes associated with porphyry ore deposits, which may represent an endmember in rapid ascent of felsic magmas from depth. In particular, we simulate the propagation of dikes launched from a deep (900 MPa), volatile-saturated magma source, testing the effects of the magma H2O/CO2 content, temperature, and mass on its ascent rate and final emplacement depth. The model predicts short ascent times (hours to days), with a large increase in viscosity at shallow depth, leading to stagnation and solidification of the dike. Higher initial water content, higher temperature and larger mass of the magma in the dike promote faster propagation and shallower arrest. Volatile loss from ascending magma remains limited until the stagnation depth, providing a potential mechanism for transfer of deep volatiles to hypabyssal blind intrusions associated with porphyry ore deposits. Our findings are applicable to the problem of silicic magma ascent through the crust more generally.

SUMMARYImplementations of Markov chain Monte Carlo (MCMC) methods need to confront two fundamental challenges: accurate representation of prior information and efficient evaluation of likelihood functions. The definition and sampling of the prior distribution can often be facilitated by standard dimensionality-reduction techniques such as Principal Component Analysis (PCA). Additionally, PCA-based decompositions can enable the implementation of accurate surrogate models, for instance, based on polynomial chaos expansion (PCE). However, intricate geological priors with sharp contrasts may demand advanced dimensionality-reduction techniques, such as deep generative models (DGMs). While suitable for prior sampling, these DGMs pose challenges for surrogate modeling. In this contribution, we present a MCMC strategy that combines the high reconstruction performance of a DGM in the form of a variational autoencoder (VAE) with the accuracy of PCA-PCE surrogate modeling. Additionally, we introduce a physics-informed PCA decomposition to improve accuracy and reduce the computational burden associated with surrogate modeling. Our methodology is exemplified in the context of Bayesian ground-penetrating radar (GPR) travel-time tomography using channelized subsurface structures, providing accurate reconstructions and significant speed-ups, particularly when the computation of the full-physics forward model is costly.

SUMMARYThis study focuses on the 3-D velocity structure and thickness of ∼7 Myr-old oceanic crust surrounding borehole 504B, located ∼235 km from the intermediate-spreading Costa Rica Rift (Panama Basin). It investigates how well seismic structure determined by 3-D tomography compares with actual lithology and, consequently, what the origin and cause might be of an amplitude anomaly, the 2A Event, that is observed in multichannel seismic data. Our P-wave model shows an ∼0.3 km-thick sediment layer of velocity between ∼1.6-1.9 km s−1 (gradient 1.0 s−1), bound at its base by a velocity step to 4.8 km s−1 at the top of oceanic crustal Layer 2. Layer 2 itself is subdivided into two main units (2A and 2B) by a vertical velocity gradient change at 4.5 km depth, with a gradient of 1.7 s−1 above (4.8-5.8 km s−1) and 0.7 s−1 below (5.8-6.5 km s−1). The base of Layer 2, in turn, is defined by a change in gradient at 5.6 km depth. Below this, Layer 3 has a velocity range of 6.5-7.5 km s−1 and a gradient of ∼0.3 s−1. Corresponding S-wave igneous layer velocities and gradients are: Layer 2A, 2.4-3.1 km s−1 and 1.0 s−1; Layer 2B, 3.1-3.7 km s−1 and 0.5 s−1; Layer 3, 3.7-4.0 km s−1 and 0.1 s−1. The 3-D tomographic models, coupled with gravity modelling, indicate that the crust is ∼6 km thick throughout the region, with a generally flat-lying Moho. Although the P- and S-wave models are smooth, their velocities and gradients are remarkably consistent with the main lithological layering subdivisions logged within 504B. Thus, using the change in velocity gradient as a proxy, Layer 2 is interpreted as ∼1.8 km thick and Layer 3 as ∼3.8 km thick, with little vertical variation throughout the 3-D volume. However, the strike of lateral gradient variation is not Costa Rica Rift-parallel, but instead follows the orientation of the present-day adjacent Ecuador Rift, suggesting a reorientation of the Costa Rica Rift spreading ridge axis. Having determined its consistency with lithological ground-truth, the resulting P-wave model is used as the basis of finite difference calculation of wave propagation to find the origin of the 2A Event. Our modelling shows that no distinct interface, or transition, is required to generate this event. Instead, it is caused by averaging of heterogeneous physical properties by the seismic wave as it propagates through Layer 2 and is scattered. Thus, we conclude that the 2A Event originates and propagates exclusively in the lower part of Layer 2A, above the mean depth to the top of the dykes of Layer 2B. From our synthetic data we conclude that using the 2A Event on seismic reflection profiles as a proxy to determine the Layer 2A/2B boundary's depth will result in an overestimate of up to several hundred metres, the degree of which being dependent on the specific velocity chosen for normal moveout correction prior to stacking.

SummaryA large non-double-couple component of a tectonic earthquake indicates that its rupture likely was complex and likely involved multiple faults. Detailed source models of such earthquakes can add to our understanding of earthquake source complexity. The 2007 Martinique earthquake in the Caribbean Sea is one of the largest recent earthquakes with a known large non-double-couple component. It was an intermediate depth intraslab earthquake within the South American plate where it is subducting beneath the Caribbean plate. We applied potency density tensor inversion (PDTI) to teleseismic P waves generated by the 2007 Martinique earthquake to model its source processes and focal mechanism distribution. We identified two focal mechanisms: a strike-slip mechanism with a north–south tension axis (T axis), and a down-dip extension (DDE) mechanism with an east–west T axis. Rupture by the DDE mechanism was predominant in the northern part of the source region and strike-slip rupture in the southern part. These two focal mechanisms had approximately parallel pressure axes (P axes) and approximately orthogonal T axes. The seismic moments released by both types of rupture were almost equal. These results indicate that the 2007 Martinique earthquake had a large non-double-couple component. We identified five sub-events with two predominant directions of rupture propagation: two strike-slip sub-events propagated to the southeast and three DDE sub-events propagated to the east. Although the directions of propagation were consistent for each focal mechanism, each sub-event appears to have occurred in isolation. For example, the rupture of one DDE sub-event propagated from the edge of the source region back towards the hypocentre. Complex ruptures that include multiple sub-events may be influenced by high pore fluid pressure associated with slab dehydration. Our results show that PDTI can produce stable estimates of complex seismic source processes and provide useful information about the sources of complex intermediate depth intraslab earthquakes for which fault geometry assumptions are difficult.

SummaryContinuous Plinian eruptions often excite atmospheric modes of ∼3.7 and ∼4.4 mHz, which are observed as harmonic oscillations of ionospheric total electron content (TEC) by global navigation satellite system (GNSS) receivers. Such TEC oscillations started shortly after the great eruption of the Hunga-Tonga Hunga-Ha'apai (HTHH) submarine volcano at ∼4:14 UT, on 2022 Jan. 15. Here I analyze GNSS data at stations within ∼4000 km from the volcano to study temporal and spatial distribution of such atmospheric modes. Strong ∼3.7 mHz TEC oscillations in near fields started shortly after the eruption onset and propagated outward with the sound speed from HTHH. Later such TEC oscillations became strong again with the amplitude peak at the distance ∼1400 km from HTHH. Such far field oscillations occurred also above New Zealand and the Solomon Islands, ∼3000 km from HTHH. Their amplitudes seem correlated with those of the 0S29 solid earth mode, suggesting that vertical surface vibrations underneath may play a role in maintaining the atmospheric mode. Onset of the far field TEC oscillations are synchronized with the local sunrise, possibly controlled by diurnal changes in the ionospheric electron density.

SummaryEarthquake source parameters can be estimated using seismological observations, but the identification of the fault responsible is often complicated by location uncertainties and the inherent ambiguity between nodal planes. Satellite InSAR can be used to observe ground deformation and model fault geometry but is limited by climate conditions (water vapour) and ground coverage (dense vegetation). In the tropics, the atmosphere is dynamic and most regions are densely vegetated, making detecting co-seismic deformation challenging. Here we perform a systematic inspection of co-seismic interferograms from Sentinel-1 SAR images, to assess their suitability for detecting co-seismic deformation in Costa Rica. Using data from the seismological network, we target seven earthquakes between 2016 and 2020 with depths ≤ 20 km and magnitudes Mw5.3 to Mw6.2 . For each event, we use the seismic parameters to compute line-of-sight displacements for ascending and descending geometries and for both nodal planes and generate 12 and 24-day co-seismic interferograms where available. We obtain interferograms with co-seismic displacement signals for three of the seven earthquakes. We invert the geodetic data to retrieve the earthquake source parameters but the lack of offshore geodetic coverage causes tradeoffs between parameters and large uncertainties. The Jacó and Golfito earthquakes likely occurred on the subduction interface and the geodetic locations were 6 − 9 km closer to the coast than previous seismic estimates. The Burica earthquake occurred on a shallow steeply dipping thrust fault in the outer forearc. For the other earthquakes, no co-seismic deformation was detected due to atmospheric noise or poor coherence. These results demonstrate the suitability of 12-day Sentinel-1 interferograms for monitoring shallow earthquakes of magnitude >Mw5.7 in Central America. This approach can be used to begin a surface deformation catalogue for the region, which will ultimately help improve the understanding of active deformation processes and improve hazard maps.