The 2021 $Mw$ 5.6 Earthquake in the Guyana Shield, South America: An Intraplate Event with No Reason to Occur Open Access

On 31 January 2021, an $Mw$ 5.6 earthquake occurred in the “Deep South” of Guyana, felt in most of the country, from Georgetown, 480 km north, to Manaus, Brazil, 650 km south. Such long distances are typical of low attenuation in intraplate settings and consistent with the magnitude–intensity relation of Quadros et al. (2019) that predicts III–IV and III modified Mercalli (MM) intensities at those distances. Severe damage to weak adobe houses and some cracks in masonry buildings (VI MM) were observed 10–20 km from the epicenter (Fig. S1, available in the supplemental material to this article). Several ground cracks were also observed at the uninhabited epicentral area, raising concernsin a country not used to earthquakes.

The January 2021 event was the largest event in the Paleoproterozoic Guyana shield. Two previous earthquakes (⁠$mb$ 4.3 in 1964 and 4.8 in 1965, 10 and 70 km west of the 2021 epicenter) suggest a local seismic zone within an otherwise aseismic region. One common explanation for intraplate earthquakes is crustal weakness zone, such as failed rifts (e.g., Schulte and Mooney, 2005; especially for large earthquakes). The Mesozoic Takutu graben is nearby (Teitz, 1991, Fig. 1) and could perhaps explain the 1965 Guyana event (Fig. 1). Errors of teleseismic epicenters in the 1960s could be as high as 30 km so that the 1965 event was classified as “possibly rift‐related” by Schulte and Mooney (2005). Or else, it could have occurred in the same fault zone of the 2021 event. However, epicentral relocation of the 1965 event (Fig. 1, Supplement Tables T1 and T2) confirmed it was farther from the 2021 epicenter and farther from the Takutu graben. Interferometric Synthetic Aperture Radar (InSAR) interferometry images and aftershock locations from a temporary local network place the 2021 epicenter 65 km south of the Takutu graben, within the Paleoproterozoic shield, hence in a “non‐rifted” continental interior. Given that the causes of intraplate seismicity are highly debated (e.g., Calais et al., 2016), a detailed study of this event contributes to that discussion and may eventually be used to test models of intraplate seismicity.

We present the results of a four‐station local network deployed in 7 March (35 days after the mainshock), and discuss source parameters such as rupture geometry, fault size, and stress drop. We also test some common models of intraplate seismicity trying to explain this rare event.

Fore‐ and aftershocks were studied with template matching at the closest permanent station BOAV, 114 km away (Fig. 1, Supplemental Text 1 and Fig. S2). Seismicity began in early December 2020, with a small peak on 25–31 December, with magnitudes up to $ML$ 3.7, felt in the local villages (Fig. S1). Foreshocks decreased during January and a few immediate foreshocks were detected preceding the mainshock of 31 January (Fig. S2b). No significant difference was found in the b‐value of foreshocks, 0.83 ± 0.08, and aftershocks, 0.86 ± 0.02 (Fig. S2a).

The mainshock magnitudes were $mb$ 5.9, $Ms$ 5.4, and $Mw$ 5.6. Table T3 lists the magnitudes given by different agencies. The regional Brazilian magnitude (Assumpção et al., 2022) was $mR$ 5.45. The regional $mR$ is equivalent to the teleseismic $mb$⁠, on average, but $mb$ tends to be higher than $mR$ for dip‐slip earthquakes (Assumpção et al., 2022), like this Guyana earthquake.

Moment tensor inversions by several international agencies (summarized in Fig. S3) consistently showed a reverse fault with southwest–northeast‐oriented P axis. A large non‐DC component, common to all inversions, ranged from 25% (Global Centroid Moment Tensor [Global CMT]) to 58% (GeoForschungsZentrum [GFZ]). The best‐fitting DC solutions of those inversions have one nodal plane striking roughly southeast–northwest and dipping southwest 45°. The hypocentral depth was not well constrained by the Global CMT inversions (nor by the teleseismic epicentral determinations), and ranged between 4 and 12 km.

A regional moment tensor was inverted using the six closest stations with good azimuthal coverage (Fig. 2). We used the ISOLA code (Zahradník and Sokos, 2025). A 76% calibrated linear vector dipole (CLVD) component was obtained, larger than found by other international agencies. Statistical analyses using a Bayesian bootstrap resampling (Fig. S4b) showed that the CLVD component is highly significant. A shallow depth can be inferred from the InSAR pattern (Fig. 3) and is confirmed by the best correlation of the waveforms at about 2 km (Fig. S4), consistent with the aftershock hypocenters determined by the four‐station local network (Fig. 4).

The large CLVD component of all teleseismic and our regional CMT inversions likely results from a complex rupture geometry (e.g., Palgunadi et al., 2020). The aftershock distribution, shown below, confirms the complex nature of the mainshock rupture.

We processed InSAR data from the Sentinel‐1 constellation on both ascending and descending tracks covering the epicentral region (track 83 with acquisitions on 31 January and 12 February and track 164 with acquisitions on 30 January and 11 February). Clear interferometric fringes highlight surface deformation (Fig. 3a). This region of Guyana is characterized by savanna‐like vegetation with little forest (except along riverbanks) and predominance of clear ground, allowing for good InSAR coherence. The area to the southwest of the magenta line (interpreted as the “surface trace”) has 14 fringes giving about 40 cm of displacement toward the satellite. Assuming ground deformation can be represented by only two directions of displacement (horizontal at azimuth N45° and vertical) the vertical motion would reach about 40 cm for 20–25 cm of notheast–southwest shortening (Fig. 3b). Only a few centimeters of subsidence to the northeast is inferred. This pattern of ground displacement indicates a predominantly reverse‐faulting mechanism with the largest surface displacement in the hanging wall right on top of the rupture centroid.

The extent of the fringes (Fig. 3a) indicates a rupture length of about 5–7 km. Figure 3a shows the epicenters (yellow and blue circles, for shallower and deeper events) of the first two days of operation of the full four‐station local network. Assuming that some of the aftershocks, recorded more than one month after the mainshock, probably occurred outside the main rupture area (e.g., Neo et al., 2021), the rupture length was about ∼5 km long.

No clearly visible trace of the rupture (such as a scarp) was found in the field, and the interferometric data does not allow us to infer precisely the rupture trace due to decoherence of the InSAR signal. However, several ground cracks were observed near the expected surface break of the fault, as exemplified in Figure 3c. These were mainly extensional cracks probably resulting from surface deformation (stretching) in the sediments and soil (only one case of a small step was observed, suggesting some component of shortening). The crack orientation (Fig. 3c) is roughly parallel to the local trend of the interferometric fringes, suggesting a relationship to the gradient of the ground displacement. It is likely that the surface trace of the mainshock rupture followed the area of highest decorrelation, as observed for larger surface rupturing earthquakes (e.g., Ross et al., 2019). A preliminary trace of the fault is tentatively shown in Figure 3a as a magenta line. The expected surface rupture is roughly along a river valley (Fig. S5) in which soft sediments may have hindered the production of a clear scarp. Finally, if the rupture did not reach the surface, it stopped very close because we do not see a clear reversal of the phase gradient near the region of high decoherence.

The vertical and horizontal displacements were tentatively modeled with a single rupture 5 km long (southeast–northwest) and 3 km wide (down‐dip), using the Okada (1985) dislocation equations (Fig. S6). A rupture plane dipping 45° to the southwest requires a slip of about 1 m to produce the observed 40 cm vertical displacement. Different dip angles and centroid positions were tried (Fig. S6a–f). The fit to the horizontal displacements is not as good as the vertical. A plane dipping 30° seems to fit slightly better the horizontal displacements. We also tried to model the surface displacement with two rupture segments. Figure S6g tried a secondary segment in the west‐northwest continuation of the main rupture, and Figure S6h tried a secondary segment with strike and dip like those of the 26 March aftershock (Fig. S8) and the minor DC from the CMT decomposition (Fig. 4). Although it is not difficult to model the vertical displacement, the horizontal displacements are harder to model.

We believe the worse fit to the horizontal displacement is mainly due to the complex rupture geometry, and the limitations of our simplistic single rectangular rupture with uniform slip. The InSAR fringe pattern suggests a rupture with orientation changing from east‐southeast–west‐northwest to south‐southeast–north‐northwest (purple line in Fig. 3a), similar to the change of the river valley seen in the topography map (Fig. S5). The full interpretation of the InSAR data requires 3D modeling with variable slip. However, the different modeling tests (Fig. S6) show that a slip of about 1 m is necessary to produce the observed 40 cm vertical displacement.

Between 5 and 7 March, four local stations were deployed (Fig. 4) in a joint effort of the Guyana Geology and Mines Commission, the Guyana Civil Defense Committee, and the University of São Paulo in Brazil. Two short‐period stations (GUY1 and GUY2) recorded at 100 samples per second, and two posthole broadband stations (3 and 4) at 200 samples per second. Station GUY1 was installed on a thin layer of sediments; stations GUY2–GUY4 had the sensor sitting directly on a granitic outcrop.

During the first four weeks of operation of the local network (7 March–4 April) about 4200 events were automatically detected and located by at least three stations, with magnitudes from $ML$ 0.3 to 4.5. The first 400 events (7–22 March) were then manually revised to improve and correct all P and S arrivals. The consistency of the arrivals was checked with a Wadati diagram (Fig. S7b). These 400 manually picked events gave a $VP/VS$ ratio of 1.706 ± 0.002. A grid search was carried out for the $VP/VS$ ratio (within the range given by the Wadati diagram) and the best half‐space $VP$⁠. The best $VP/VS$ ratio was 1.705, and the half‐space $VP$ was 5.89 km/s, which produced the lowest average root mean square (rms) residual of 0.0088 s. For comparison, the average rms residual of the automatically located events was 0.02 s.

Figure 4 shows the aftershock distribution of the selected best events (P and S at all four stations, epicenter uncertainty < 60 m, and depth uncertainty < 250 m). A clear east‐southeast–west‐northwest trend of the hypocenters can be seen, consistent with the major DC component of the CMT inversion (red focal mechanism plot in Fig. 4). Besides the main ∼5‐km‐long rupture, identified by the south‐southwest‐dipping group of events in the center (inside the rectangle), other small clusters are seen to the northwest and to the southeast of the main rupture. These other clusters seem to define other smaller faults with different orientations, which could be attributed to triggered aftershocks from Coulomb stress changes around the mainshock.

For example, the northwest aftershock cluster has a different orientation steeply dipping to the north‐northeast. The largest aftershock (26 March 2021, $mb$ 4.6 and $Mw$ 4.1) occurred near this northern cluster (Fig. 3, Fig. S8). The focal mechanism determined with the FMNEAR code (Delouis, 2014) by waveform modeling of P and S waves, in the 0.4–1.0 Hz frequency range, showed a dip‐slip event with a small strike‐slip component on a fault‐plane striking east–west and steeply dipping to the north (Fig. 4, Fig. S8a), consistent with the impulsive P‐wave polarities of the four local stations and BOAV. This 26 March event had one immediate foreshock and three aftershocks for which relative location, using correlated P and S picks, showed a fault‐plane striking 285° and dipping 69° to the north (Fig. S8b) consistent with the FMNEAR focal mechanism solution and all five P‐wave polarities. Similarly, the southeast cluster suggests a secondary west‐northwest–east‐northeast‐oriented fault dipping to the north (Fig. 4).

It should be noted that the aftershocks in the central part (events enclosed in the large rectangle in the map of Fig. 4) lie along an apparent dip shallower than the 41° dip of the major DC component (seen in the projection plot labeled “Main rupture”). One possibility to explain this inconsistency is that most aftershocks could occur inside the hanging wall and not on the plane of the main rupture. Similarly, some aftershocks to the northeast of the main rupture (events in the right side of the “Main rupture” projection) could have occurred in the footwall, also off the main rupture.

The aftershocks near the main rupture (central rectangle in Fig. 4) do not seem to align along a single planar plane. Although some of the events may be off the main rupture, the scatter of the hypocenters suggests a complex geometry of the mainshock rupture. This complex rupture geometry is likely the best explanation for the high CLVD component of all CMT solutions.

The shallow nature of the rupture and the unexpectedly large observed ground displacement make this event of special interest to studies of intraplate seismicity. Assuming the rupture size is about 5 km long (which is consistent with the $Mw$/rupture‐length empirical relation of Wells and Coppersmith, 1994 and Ciardelli and Assumpção, 2019) and 2.5–3 km wide down‐dip (Fig. 4), a 0.8–1.0 m average slip is required to produce the ∼40 cm vertical ground displacement (see Fig. S6 and Text S2). Surprisingly, no obvious trace of the fault at the surface was recognized. An average slip of 1 m is larger than expected for an $Mw$ 5.6 event and is more common for M ∼7 earthquakes but is still near the observed range for reverse‐faulting earthquakes (Wells and Coppersmith, 1994).

A rupture area of $5×2.75km2$ with 1 m average displacement gives a moment magnitude of $Mw$ 5.7, consistent with the moment tensor inversions and the empirical $Area×M0$ relation for stable continental region earthquakes (Leonard, 2014). Relating the rupture area with the estimated slip and $Mw$ gives a stress drop in the range 5–12 MPa, larger than the average for shallow events in midplate South America (Ciardelli and Assumpção, 2019). Text 2 of the supplemental material discusses these estimates of stress drop, their uncertainties, and implications.

The shallow nature of the 2021 $Mw$ 5.6 Guyana earthquake seems to be typical of other midplate earthquakes in South America, which tend to occur mostly in the upper 5 km of the crust (Assumpção et al., 2016). In eastern and central Brazil, moderate‐size earthquakes have occurred with well‐determined shallow depths, such as the Mara Rosa 2010 earthquake (⁠$Mw$ 4.4 at 1.4 km depth; Barros et al., 2015) and the Itacarambi 2007 event in the São Francisco craton (⁠$Mw$ 4.7, 0.7 km depth; Chimpliganond et al., 2010), neither of which was anthropogenic. Similar shallow depths have been observed in the cratonic region of southwest Australia for magnitudes in the $Mw$ 4.4–6.0 range (e.g., Dawson et al., 2008).

Despite the shallow depth and the relatively flat topography of the epicentral region, no prior scarp was recognized in the valley so no indication on the surface of repeated strong earthquakes in the recent geologic record.

The regional stress field in the Amazon craton is not well defined yet due to the few available focal mechanisms (Assumpção et al., 2016). In the southern part of the craton, strike‐slip mechanisms predominate with northwest–southeast‐oriented compression. In the central part of the craton (Amazon basin and Guyana shield) reverse faults predominate with an average P axis in the northwest–southeast to north–south direction. The mechanisms of the 2021 Guyana mainshock and its large 26 March aftershock have northeast–southwest‐trending P axes, which are not entirely inconsistent with the other nearby earthquakes in the central part of the craton, but may hint to a slight change in the orientation of the compressional stresses in that part of the Guyana shield.

Explanations for intraplate seismicity are still highly debated, and several different models have been proposed (e.g., Mazzotti, 2007). Most models include either a crustal weak zone or a mechanism of stress concentration. The most studied candidates for weak zones are Mesozoic failed rifts and passive continental margins revised by Schulte and Mooney (2005). In fact, given the epicentral uncertainty, the 1964 and 1965 events had been classified by Schulte and Mooney as possibly related to the Takutu graben. As discussed earlier, the 2021 Guyana earthquake, as well as the earlier 1965 event, occurred more than 60 km from the Takutu graben (Fig. 1). In addition, all major faults in the epicentral Paleoproterozoic terrains trend southwest–northeast, opposite to the northwest–southeast orientation of the 2021 rupture. So, the Guyana earthquakes cannot be attributed to any weak zone, or reactivation of old fault zones. It seems to have occurred in a minor secondary fault not shown in regional geological maps.

Stress concentration due to lithospheric thinning has been proposed to contribute to earthquakes in Brazil (Assumpção et al., 2004), and other regions (Lebedev et al., 2023). Here, we compare the seismicity in midplate South America, given by the catalog of the Brazilian Seismic Network (Bianchi et al., 2018), updated to 2023, with two models of lithospheric S‐wave velocities at 110 km depth: model SAAM23 (Ciardelli et al., 2022) and model SACI24 (Melo et al., 2024). S‐wave anomalies at this depth are a proxy for lithospheric thinning. As shown in Figure 5a,b, both models indicate that in most of the study area, there are more epicenters in areas of low‐$VS$ anomaly (lower than +1%) than should be expected if the epicentral distribution had no relation with lithospheric thinning. However, in the Guyana epicentral area, there is no low‐$VS$ anomaly (shown by the stars along the axes of the histograms in Fig. 5a,b). This means that lithospheric thinning cannot explain the Guyana seismicity.

Stress concentration in the upper crust, due to flexural deformation caused by lithospheric load, explains seismicity in central Brazil and a few other areas (Assumpção and Sacek, 2013). Lithospheric load can be inferred from free‐air gravity anomalies. We use two different anomaly maps, one by Sá (2004), based mainly on terrestrial measurements, and the gravity field and ocean circulation explorer (GOCE) satellite map. Figure 5c,d shows that, in mid‐plate South America, areas with positive gravity anomalies have more earthquakes than would be expected if seismicity were uniformly distributed. However, at the epicenter of the Guyana earthquakes, there is no positive gravity anomaly that could indicate a lithospheric load causing flexural stresses in the upper crust.

In the isolated, scarcely populated epicentral area of the “Guyana Deep South”, there is no mining activity, nor any hydroelectric reservoir that could possibly induce earthquakes. In addition, stresses from icesheet unloading from the last glacial maximum are thought to contribute to seismicity in high‐latitude regions, which is not the case for Guyana. Therefore, none of the most common models for intraplate seismicity seem to be applicable to the Guyana earthquakes. This has implications for seismic hazard studies because it reinforces the challenge of using geological and geophysical data to define seismic zones. On the other hand, the 1964, 1965, and 2021 Guyana earthquakes all seem to belong to a large seismic zone that stretches as farther south as latitudes of 5° S, near the Amazon River (Fig. 5). This means that identifying clusters of epicenters is still an important criterion in seismic hazard studies, even when the seismogenic mechanisms are not known yet.

The tomographic maps of Figure 5 seem to indicate a concentration of epicenters in the edges of cratonic deep keels, such as suggested by Mooney et al. (2012) on a global scale for magnitudes larger than $Mw$ 4.5. This motivates more detailed quantitative studies in northern South America in the future.

The 31 January 2021 $Mw$ 5.6 Guyana earthquake, the largest event recorded in the Guyana shield, occurred on a shallow reverse fault with a centroid about 2 km deep. Although several extensional cracks have been observed in the epicentral area, no clear sign of a rupture scarp was identified. The reverse‐faulting mechanism is typical of the central part of the Amazon craton, but the 5–12 MPa stress drop is higher than other shallow earthquakes in mid‐plate South America.

Aftershock locations suggest a complex rupture with a main northwest–southeast striking, southwest‐dipping plane, possibly steepening toward the northwest and southeast sides. This complex rupture explains the large CLVD component of the regional (Fig. 2) and global moment tensor inversions (Fig. S2).

No clear mechanism was found yet to explain this mid‐plate activity. The regional geology (exposed Paleoproterozoic rocks), little sedimentary cover, and flat topography do not indicate any crustal weakness: no indication of fault reactivation was found as all the old basement faults strike southwest–northeast, opposite to the rupture orientation. This makes the Guyana earthquake another example of intraplate seismicity not associated with a weakness zone (failed rift or paleosuture). In addition, other models of stress concentration used to explain intraplate seismicity in other parts of South America (lithospheric thinning and flexural stresses from lithospheric loads) cannot be applied to the 2021 Guyana earthquake.

Parametric data of the mainshock and largest aftershocks were taken from the International Seismological Centre (ISC) Bulletin (International Seismological Centre [ISC], 2024) and the U.S. Geological Survey, and waveform data from the Incorporated Research Institutions for Seismology (IRIS) database. Data from the Brazilian Seismographic Network (RSBR) is open: http://www.sismo.iag.usp.br/rq/. Data from the four‐station local network are available upon request. Satellite Interferometric Synthetic Aperture Radar (InSAR) data were obtained from the European Space Agency, now available at https://sentinels.copernicus.eu/. All websites were last accessed in June 2025. The supplemental material has more detailed data.

The authors acknowledge that there are no conflicts of interest recorded.

The authors thank the Guyana Consulate in Boa Vista, Roraima, for efficiently issuing visas and permits during the COVID‐19 pandemic. The Civil Defense Commission, Guyana, provided important support for the field work. Ricky Singh and Ibrahim Abudu (GGMC) are thanked for their work on the field survey of the ground effects. Bruno Collaço and José Roberto Barbosa analyzed regional records for the Brazilian Seismic Bulletin. Francis Mandakin (Katoonarib Council), is thanked for essential help in the field work. M. A. was supported by CNPq Grant 30.1284/2017‐2, and FD by Petrobras Grant 2017/00159‐0. E. C. and R. J. acknowledge support from the Institut Universitaire de France and the European Research Council under the European Union’s Horizon 2020 research and innovation program (Grant 101125232). The authors thank Bertrand Delouis for the FMNEAR code. Emilia Brasilio and Bruno Collaço helped implementing the Okada (1985) routines provided by Scott Henderson. John Adams, John Ebel, and two other anonymous reviewers are thanked for helping improve this article.