On the Origin of Orphan Tremors and Intraplate Seismicity in Western Africa

On September 5–7, 2018, a series of tremors were reported in Nigeria’s capital city, Abuja. These events followed a growing list of tremors felt in the stable intraplate region, where earthquakes are not expected. Here, we review available seismological, geological, and geodetic data that may shed light on the origin of these tremors. First, we investigate the seismic records for parent location of the orphan tremors using a technique suitable when a single-seismic station is available such as the Western Africa region, which has a sparse seismic network. We find no evidence of the reported tremors within the seismic record of Western Africa. Next, we consider the possibility of a local amplification of earthquakes from regional tectonics, reactivation of local basement fractures by far-field tectonic stresses, post-rift crustal relaxation, landward continuation of oceanic fracture zones, or induced earthquakes triggered by groundwater extraction. Our assessments pose important implications for understanding Western Africa’s intraplate seismicity and its potential connection to tectonic inheritance, active regional tectonics, and anthropogenic stress perturbation.


INTRODUCTION
On September 5-7, 2018, a series of low-magnitude tremors hit Nigeria's capital city in Mpape, Abuja (Government Report, 2018). These events followed a series of earthquakes felt in the region since 1933, a stable, intraplate setting, otherwise not being earthquake-prone (Akpan and Yakubu, 2010;Akpan et al., 2014;Tsalha et al., 2015). This seismicity highlights the issue of seismic hazards in the Western Africa region and the need for an improved seismic monitoring network (Alaneme and Okotete, 2018;Afegbua et al., 2019). Nigeria is located in the southern part of the Neoproterozoic Trans-Saharan Mobile Belt, separating the Archean West African Craton, Congo Craton, and the Archean-Proterozoic Sahara Metacraton. Within this mobile belt, a large continent-scale system of elongate rift basins (aulacogens) developed during the Cretaceous, among which is located the Benue Trough on whose flank the Abuja city is located ( Figure 1A). The extent of areas that historical seismicity felt in Nigeria encompasses regions within the failed rifts and areas of the exposed basement on the flanks of the rift basins ( Figure 1A). Note that in the context of this paper, the term "tremor" is used in its general meaning with no particular emphasis on the processes that are attributed to a tremor activity, such as fluid migration and slow slip events, observed in some volcanic and tectonic settings (Rogers and Dragert, 2003;McCausland et al., 2005;Nadeau and Dolenc, 2005;Wech et al., 2012;Montgomery-Brown et al., 2013).
Compared to historical events, the series of shakings reported on September 2018 is located within the northern edge of Abuja. Because the capital city of Abuja is highly populated and relatively affluent, the considerable shaking (which was felt by most of the residents) was reported by the local news agencies and picked up by increased social media activity. Although the local population felt these events, very little observational evidence exists for the origin of the shaking and its connection to local geology, or regional tectonics. Here, we present the available geological and geophysical constraints that may offer clues about the shaking and examine the available seismic, geologic, and geodetic measurements and discuss several viable hypotheses regarding the origin of the felt shaking connected to the local or regional tectonics and anthropogenic activities. In particular, we evaluate the hypotheses that these events are due to 1) local amplification of earthquakes from regional tectonics, 2) distant teleseismic events large enough to be felt in Western Africa, 3) reactivation of local basement fractures by far-field tectonic stresses, post-rift crustal relaxation, or landward continuation of oceanic fracture zones, or 4) local anthropogenic activities such as groundwater extraction.
A comprehensive search of the global earthquake catalog in the decade leading up to September 2018 does not turn up any events located within Nigeria, although we point out that in the Western African region, these catalogs are incomplete below moment magnitude 4 (Kadiri and Kijko, 2021). In the region, most of the events are clustered around oceanic transform faults, with only a few inland earthquakes. Similarly, all earthquakes in 2018 occurred along the oceanic transform fault, except for a single event on February 2019, located in the country of Guinea about 10°northwest of Nigeria ( Figure 1). Therefore, we examine the available seismic data from the Nigerian seismic network (Afegbua et al., 2011), supplementing them with seismic stations from the Western African region where long-term real-time monitoring is available. We conduct a rigorous search of the seismic record using a single-station detector method to examine the first two hypotheses, which investigate the possibility of the shaking being related to distant earthquakes and their ability to reactivate preexisting fault systems that may be preferentially primed for failure (Han et al., 2017;Neves et al., 2018). We then complement the seismic investigation with other geological and geophysical datasets. Surface geology and basement structure are FIGURE 1 | (A) West African seismicity and tectonics, showing the epicenters (circles) for earthquakes in the past decade with M3.0 and greater, ocean fracture zones (red lines), and the seismic network (triangles) available for seismic detection. The 2018 earthquake events (red circles) are outside the Abuja area (red polygon north-east of TORO). OTFZ oceanic transform fault zone (oceanic fracture zone). Tectonic structures are from Pérez-Díaz and Eagles (2014). The Mesosoic failed rifts shown in yellow lines define the West African Rift System. (B) Final selection of five stations (green triangle) used to search for an earthquake with a hypothetical epicenter (red star) located in Nigeria (NG). Gray stations do not pass the data quality selection criteria. Note that although TAM (S5) has good quality data, it is not used for single station detection since it only has one horizontal channel operating during the time of interest. investigated using high-resolution satellite images and aeromagnetic data, and the observed patterns compared to the inferred prevailing stress field and the general trends dictated by the opening and closure of the failed rifts located close to the Abuja capital city. A final analysis explores the anthropogenic controls on earthquake nucleation and whether these may be related to the hydrological cycle within Nigeria. We measure surface deformation using Interferometric Synthetic Aperture Radar (InSAR) and correlate the observed patterns with the ongoing anthropogenic activities such as groundwater pumping and the addition of dams in the study area. Overall, the assessments presented in this study pose important implications for understanding Western Africa's intraplate seismicity and its connection to regional tectonics and local geology.

Single-Station Seismology
Ideally, the detection of a low-magnitude local or regional event would require a proximal high-quality seismic array. Western Africa, on the other hand, has a sparse distribution of real-time seismic stations. Even recently deployed small-aperture seismic arrays located in Western Africa, e.g., Nigeria (Afegbua et al., 2011) and Ghana (Ahulu and Danuor, 2015), do not currently archive data on global waveform databases and only a few are real-time. Therefore, we use a single-station location method that is based on polarization analysis using eigen-decomposition of ground displacement (Vidale, 1986;Park et al., 1987;Bai and Kennett, 2000;Simons et al., 2009) for event detection, seismic phase identification (Earle, 1999), and source localization (Böse et al., 2017). Alternative methods such as beamforming analysis, which is widely applied to pin-pointing earthquake source direction (Rost, 2002;Rost et al., 2006;Nakata et al., 2019), are challenging to apply because of the sparsity of high-quality, continuously recording, and small-aperture seismic arrays present at regional distances to the location of the largest reported shaking.
The single-station location technique is based on the same principle for which polarization analysis is most commonly used, i.e., orienting seismometers on the bottom of the seafloor (Stachnik et al., 2012;Zha et al., 2013;Doran and Laske, 2017;Scholz et al., 2017) or identifying misoriented horizontal channels (Ojo et al., 2019). Used in this mode, the arrival direction (azimuth) and the distance of a seismic event can be inferred from a three-channel seismogram recorded on a single station on the African continent. The geographical orientation of the seismic channels is then determined, given a known arrival direction of a particular seismic phase (often the compression and Rayleigh waves). In the earthquake-location mode, however, the idea is reversed. We assume the channel orientations are known (or can be corrected for); then, in principle, we fix the source orientation by determining the azimuth of the incoming waves. The key idea is to identify the polarization of particular seismic phases which are parallel to the wave propagation (e.g., compressional and Rayleigh wave phases) and use this to determine the wave-arrival azimuth. A differential time between the first arriving compressional wave and other phases (e.g., S, Love, and Rayleigh), in conjunction with a given earth model, prescribes the epicentral distance and completes the process, providing coordinates of the earthquake.
The adaptation of polarization analysis for single-station event location (Magotra et al., 1987;Magotra et al., 1989;Frohlich and Pulliam, 1999;Agius and Galea, 2011;Böse et al., 2017) increases the likelihood of detecting earthquakes with smaller magnitudes. In our adaptation, we incorporate benefits from the coherence (Vidale, 1986) and covariance technique (Jurkevics, 1988;Schulte-Pelkum et al., 2004;Park and Ishii, 2018) and include multimode identification (Bai and Kennett, 2000;Bai and Kennett, 2001). Our idea bears some similarity to the technique proposed for locating Mars-quakes (Böse et al., 2017), or the identification of the source location of ambient seismic noise or micro-tremors (Koper and Hawley, 2010;Zha et al., 2013). We describe a successful application for locating the epicenter of moderate-to-large teleseismic earthquakes. We then describe the application to a few high-quality stations recording during the 3 days with the largest shaking, i.e., stations on the GEOSCOPE network (Roult et al., 2010), the Global Seismic Network (GSN) (Lay et al., 2002), and the Nigerian Seismic Network.
While single-station event location is sufficient for distinguishing between a teleseismic and a local/regional event, when more than one station is available, the uncertainty or bias (error) from a single-station detection can be improved by performing a second-stage association analysis (Ringdal and Husebye, 1982;Shearer, 1994;Ekström, 2006). This allows us to improve the confidence estimate for the location by comparing the estimated source locations from multiple stations when they are available. Suppose all the single-station locations agree (within uncertainty) that the hypothesized source location (Abuja, Nigeria) is not the epicenter of the detection, in that case, we can rule out this null hypothesis and explore the alternative hypothesis that the shaking is from a different source location (i.e., of teleseismic origin). We require only two stations to confirm detection for a particular source location. Application to the 360 polarization records derived from each of the five closest stations to Nigeria results in a space-time probability of a possible event detection (Böse et al., 2017) (see Supplementary Appendix for a detailed discussion of technique).

Structural Mapping From Aeromagnetic Data Analysis and Satellite Images
The largest shaking may have resulted from seismogenic slip on local faults or due to local amplification. To explore how local geology in and around Abuja may have contributed to significant shaking intensity, following standard practice, we analyze the subsurface structure and the spatial distribution of the sedimentary cover. The Abuja area is located in low magnetic latitude. Thus, we first perform a reduction of the magnetic data to the magnetic equator (RTE, e.g., Li and Oldenburg, 2001) using an IGRF 2005 model. This transformation corrects for the skewness of the magnetic anomalies due to the oblique angles of magnetization at large distances from the magnetic pole.
We then upward-continue the RTE-corrected aeromagnetic grid to remove noise and apply a vertical derivative filter to better resolve the gradients that correspond to the geological structure. Further, we estimate the distribution of depths to the top of the magnetic sources (i.e., crystalline basement) using the Source Parameter Imaging (SPI) method (Smith et al., 2002;Smith and Salem, 2005). Although this technique has an accuracy of ±20%, it can reliably show the relative spatial pattern of the depth of the burial of the magnetized basement (e.g., Kolawole et al., 2018). Volcanic rocks are absent in the Abuja area; thus, the Precambrian metamorphic rocks and granitic intrusives define the primary magnetic basement. In order to identify the dominant structural trends associated with basement deformation (and possibly faulting) in the area, we manually digitize the linear traces defined by edges in the filtered (vertical derivative) aeromagnetic grid. Further, we map fracture systems in basement exposures in satellite images (Google Earth©) and compare with the aeromagnetic fabric trend to better constrain the dominant basement structures that may correspond to faulting.

InSAR Deformation Field
To map the surface deformation, we use 80 SAR images acquired in the ascending orbit of the Sentinel-1A/B C-band satellites between 2018/01/19 and 2020/09/05. We performed an advanced multitemporal SAR interferometric analysis to retrieve rates and time series of surface deformation over the study area. The analysis began with co-registering SLC images to a reference image, which includes a standard matching algorithm using a digital elevation model (DEM), precise orbital parameters, and amplitude images (Sansosti et al., 2006). For the Sentinel-1A/B datasets, the step mentioned above is followed by an enhanced spectral diversity (ESD) approach (Yagüe-Martínez et al., 2016;Shirzaei and Bürgmann, 2017). Using this dataset, we generate a set of high-quality interferograms, considering only those with short perpendicular and temporal baselines. We also apply a multi-looking operator of 32 and 6 pixels in range and azimuth to obtain a ground resolution cell of ∼75 m × 75 m. To calculate and remove the effect of topographic phase and flat earth correction (Franceschetti and Lanari, 1999), we used a 1-arcsecond (∼30 m) Shuttle Radar Topography Mission DEM (Farr et al., 2007) and precise satellite orbital information. To identify the elite (i.e., less noisy) pixels, we only consider pixels with an average temporal coherence larger than 0.65. In analogy to Global Navigation Satellite System (GNSS) stations, these elite pixels form a network of SAR observation points, for which we estimate the time series and rate of line-of-sight (LOS) surface deformation. To retrieve the absolute (unwrapped) phase values for sparsely distributed elite pixels, we applied a minimum cost flow (MCF) algorithm. Although the precise orbits are used, a few interferograms were still affected by a ramp-like signal, which was corrected by fitting a second order polynomial to their unwrapped phase (Shirzaei and Walter, 2011). We further applied several wavelet-based filters to correct for effects of spatially uncorrelated topography error and topography correlated atmospheric delay (Shirzaei and Bürgmann, 2012). Subsequently, we applied a re-weighted least square approach iteratively to invert the corrected measurement of the unwrapped phase at each elite pixel and solve the time series of the surface deformation. We further reduce the effect of residual atmospheric errors by applying a high pass filter based on continuous wavelet transform to the time series of surface deformation at each elite pixel. Finally, we estimate the long-term LOS deformation rates as the best-fitting line slope to the time series of surface deformation at each elite pixel.

Earthquake Detection with a Sparse Network
We use the following station-selection criteria: 1) waveform data recorded during tremor activity archived and retrievable from a global waveform database, e.g., IRIS, 2) good back-azimuth coverage relative to Nigeria, 3) availability of three-channel seismograms necessary for polarization analysis, 4) high signal to noise ratio, and 5) correct orientations for the horizontal channels or new analysis reflecting the proper orientations (e.g., Ojo et al., 2019). Only five stations out of the ten initially selected pass our selection criteria, and their locations provide relatively good azimuthal coverage around Nigeria (see Figure 1B and Tables 1, 2). We focus on the 3 days from September 5-7, 2018, when anecdotal reports agree on the strongest shaking in Nigeria, 1) three reports on September 5: ∼13:30, 16:30, and 19:00 UTC, and 2) two on September 6 and 7: 01:30 and 05:30 UTC. We scan hour-long records of three-component seismograms at each station, allowing for a half-hour to extend detection duration. Our dataset results in 360 polarization records, representing 72-hourlong records (for the 3 days) at each of the five stations. A data quality check reveals that two of the five stations are particularly useful for seismic detection (i.e., G.TAM and G.DBIC). In particular, stations NJ.TORO, which the closest station to the main interest area, Abuja, is very noisy and is impractical for detection on its own (see Figures 1A, 2). This analysis further highlights the challenge of earthquake location and the value of a single-station event detector.

Null Detection and Other Coincidences
We analyze the continuous data stream in hour-long sections during the 3 days of the largest reported shaking. In this 3-day scan, we detect a notable event in the hour following 18:00 UTC on September 5th with good SNR on most stations ( Figure 2). Our attempt to associate the other tremor sequences to either a previously undetected global or regional event proved inconclusive (Supplementary Figures S5A-D). For the September 5 event, the stations detect P and S arrivals within only slight timing variations. We cross-reference our detection with the USGS earthquake catalog. The arrival time of the seismic phases is correctly matched with a source location that is consistent with the M6.6 Hokkaido Eastern Iburi, Japanese earthquake (Figures 3,  4B). We note that this event is neither local nor regional (<30°) but clearly of teleseismic origin (compare Figures 4A,B with Figure 3). Despite the large station-earthquake distances (80-120°), we detect clear body-wave and surface-wave phases on the high-quality stations located around Nigeria. Even for the station at the largest epicentral distance (GT-DBIC), strong PP and SS phases are visibly detected and identified by our phase detector ( Figure 3 and Table 2). The closest stations to the Japanese event are able to pick up at least one of the two direct P and S arrivals. For example, station II-RAYN demonstrates that our method provides effective picks for four of the incoming phases (P, PP, S, SS, L, and R: Figure 3). At stations MN-WDD and IU-FURI located at an epicentral distance less than 90°, we precisely identify all the main P and S arrivals, except for the reflected SS phase. At station II-MBAR, which is even further away (>90°) from the source location, our method is still able to detect the diffracted phases (Pdiff and Sdiff). However, at station GT-DBIC which is located ∼120°from Japan, we are unable to detect the first-arriving diffracted Pdiff waves and only able to detect the diffracted Sdiff phase (with a slight time shift compared to the travel time predictions using the ak135 Earth model) and the reflected PP and SS phase. Association of the single-station detections shows that three of the five stations (WDD, FURI, and MBAR) provide a location match with an error <5°, while the other two stations (DBIC and RAYN) agree on the epicentral distance with some bias in the azimuth direction ( Figure 5). This is clear evidence that the seismic event detected at all our stations is coincident in space and time with the teleseismic Japanese event reported by the global catalog and felt in Nigeria on September 5th.
Our search for a source-origin in Western Africa was inconclusive, which may suggest a low magnitude event that may have been attenuated before being detected by the regional seismic network. We emphasize that of the five events reported in anecdotal records; we are able to associate a single event to a teleseismic earthquake originating in Hokkaido, Japan, on September 5, 2018 (18:07 UTC). This event, from far-away Japan, triggered multiple landslides (∼6,000) that lasted for a few minutes (Yamagishi and Yamazaki, 2018;Kameda et al., 2019;Shao et al., 2019;Wang et al., 2019). None of the other tremors in the sequence identified by anecdotal records are conclusively associated with either a local or teleseismic event.
We explore briefly what explanation, if any, exists between the felt shaking and this landslide event. Alternatively, we evaluate how the local geology around Abuja could host the reactivation of existing faults. We also consider the possibility of induced earthquakes triggered by groundwater extraction.

Geological Framework from Aeromagnetic and Satellite Data
We integrate available surface and subsurface geological and geophysical datasets to constrain the first-order geological features in the Abuja area that could have localized or amplified ground shaking. First, we delineate satellite-scale fracture systems in basement outcrops using Google Earth 1 | Stations used to discriminate local and teleseismic events originating from Western Africa. Station G.TAM has only one horizontal channel, station II.ASCN is an island station, and station PM.PFVI is located in Portugal, but with very little propagation in the Atlantic Ocean. We include II.RAYN, which is a very quiet station, to broaden the back-azimuth coverage, even though it is technically in the Middle East. The bold indicates the five high-quality stations with very clear event detections (compare Figure 3 with Figure 2). a Proper orientation for the North Channel is retrieved from Ojo et al. (2019). H 0 : Hypothesis that event is regional with epicenter in Nigeria. H 1 : Hypothesis that event is teleseismic with epicenter in Japan.
[Δ*, Φ*]: Observed distance and azimuth (see Figure 5) [ΔH, ΦH]: Calculated distance and azimuth assuming Hypothesis (compare Figure 5 with Figure 4) ϵ: Normalized error between observation and hypothesis (low error confirms teleseismic hypotheses, H 1 ) Frontiers in Earth Science | www.frontiersin.org September 2021 | Volume 9 | Article 716630 FIGURE 2 | Seismograms at stations in or close to Nigeria using an 80 min time window. Traces are amplitude normalized. Stations used for further analysis, marked in black, are those that pass the data selection criteria (green triangles in Figure 1B). In each trace, the first P (blue) and S (red) arrivals of the M6.6 Hokkaido event are shown with the station code to the top.
Orphan Tremors in Western Africa FIGURE 3 | Phase detection using polarization analysis applied to data recorded at the five high-quality stations during the Hokkaido event of September 5, 2018. The waveforms are for an 80 min-long time window, where time 0:00 is 20 min before the P arrival at each station following the event origin. The detected arrival of P and S waves is marked with blue and red dashed lines, respectively, while the first arriving linearly polarized surface wave is shaded in gray. The predicted arrival time calculated using the earthquake parameters and an earth model is indicated with dots and labels indicating the phase.
Frontiers in Earth Science | www.frontiersin.org September 2021 | Volume 9 | Article 716630 satellite images at a spatial resolution of 5 m (e.g., Kolawole et al., 2019). Then, we utilize high-resolution aeromagnetic data first to map sub surface structural fabrics in order to obtain information on the potential trends of fault lineaments in the crystalline basement and then to model the distribution and thickness of the sedimentary over-burden (e.g., Grauch and Hudson, 2007;Kolawole et al., 2018). The aeromagnetic data used in this study were acquired country-wide between 2005 and 2007 by the Nigerian Geological Survey Agency (NGSA), with 500 m line spacing and 80 m mean terrain clearance. The aeromagnetic data were provided as a grid of 100 m cell size.

Delineated Geological Structures in the Abuja Area
The satellite-scale fracture systems in the basement outcrops, mapped in the satellite images of the Abuja area ( Figures 6A,B), show dominant trends of 049°± 3.6°(NE-SW) and 136°± 3.7°(NW-SE) ( Figure 6C). Also, the larger-scale basement fabrics, delineated in filtered aeromagnetic data ( Figure 6D), show dominant trends of 048°± 6.4°(NE-SW), 124°± 6.1°(NW-SE), and 086°± 5.4°(E-W) ( Figure 6E). Further, our depth-to-basement map, generated from the SPI transform of the aeromagnetic data ( Figures 6F,G), shows that most of the Abuja area and mainly the areas where the tremor was felt are located within a small sedimentary basin. The thickness of this sedimentary cover ranges between ∼200 and 800 m ( Figure 6G). Although the three prominent trends (NE, NW, and E-W) can be observed in the northern half of the area where basement outcrops dominate, the NE-trending lineaments appear to dominate the southern part where the basin is located ( Figure 6D). Figure 7 shows the rate of LOS displacement field and selected time series at the sites of rapid subsidence. The negative values (cool colors) correspond with movement away from satellite, hereafter, subsidence. The map is characterized by widespread subsidence up to 35 mm/yr. Roughly the zone of subsidence is bounded by the administrative divides, suggesting an anthropogenic drive. A rapidly declining trend characterizes the selected time series of LOS displacement at sites (b), (c), and (d) until 2020. The slightly rising time series ( Figure 7E) also follows a similar pattern, whose rising trend is interrupted by 2020. This behavior change might be attributed to the COVID-19 global pandemic, which has caused a reduction of economic activities worldwide.

Origin of Tremor: Perspective, Interpretation, and Outstanding Questions
In our study, we were able to associate only one of the tremor sequences to an event in the global catalog. Our attempt to associate the other events in the tremor sequences to either a previously undetected global or regional event proved inconclusive. This lack of corroborating seismic evidence confirming anecdotal reports could be due to either of the following: 1) the inaccuracy of the timing or the number of the reports; 2) localized low-magnitude shaking not detectable by the seismic station closest to Abuja, due to the quality of the data recovered from that seismic station (i.e., TORO). While this is an unsatisfactory conclusion, it emphasizes the need for high-quality FIGURE 4 | Station locations relative to two hypothetical epicenters (star). (A) For a local or regional event originating from Nigeria or within Western Africa, all stations are expected to be within 30°of the event location (red star). (B) For the Hokkaido event, distance is 60-120°and azimuth 270-320°(compare Figure 5).
Frontiers in Earth Science | www.frontiersin.org September 2021 | Volume 9 | Article 716630 data in regions that are often ignored due to the assumption of lack of seismicity. We also demonstrated that one of the reported tremor activities in Abuja, Nigeria, and recorded by the Western African stations coincides with the Hokkaido Eastern Iburi, Japanese earthquake of September 5, 2018. It is very puzzling that this moderate-sized earthquake (∼M6.7) generated enough seismic energy to be felt in the capital city of Nigeria, considering that Japan is located at considerable teleseismic distances (∼112°).
Here, we discuss a few plausible scenarios that may offer more insights into these results taken together.

Hypothesis I: A Unique Teleseismic Event
Moderate teleseismic events are unlikely to cause significant shaking at such a considerable distance. However, we are persuaded to explore the possibility that there is a direct connection to the M6.6 Hokkaido Eastern Iburi, Japanese earthquake on September 5, 2018, given a large number of coseismic landslides ) generated significant long-period surface wave energy that propagated to large distances without much dissipation (Allstadt, 2013). Landslides can be effectively modeled by single forces, where the orientation of the force exerted on the ground is in the direction opposite to landslide acceleration, generating maximum surface wave amplitudes along the trajectory of mass loss (Kawakatsu, 1989;Ekström et al., 2003;Tsai and Ekström, 2007). For shallow landslides with estimated mass volumes and a simple trajectory, it should be straightforward to model the effective forces, using synthetic seismograms (Ekström and Stark, 2013). However, the Hokkaido Iburi Earthquake triggered ∼6,000 landslides with complex trajectories (Figure 8), distributed over a large area making this a difficult task. Despite this complexity, we observe that the aggregate spatial distribution of the landslides can explain surface-wave radiation in the direction of Nigeria (310°, which is in the NW/SE direction from Japan) (Figure 8).
A comprehensive inventory of landslides Zhang et al., 2019) document important characteristics that fit this pattern: a large concentration of 65% of the landslides (21 per square-km) in an elliptical area of 173 squared-km with major axis oriented NNW/SEE at 327.7°(within 18°of the direction of Abuja). We expect that a considerable portion of the radiated surface-wave energy is connected to the coseismic landslides FIGURE 5 | Single-station location of Hokkaido, Japan, in the event reference frame (compare with station-event geometry in Figure 1B). All stations identify the correct epicentral distance and azimuth range (except for DBIC and RAYN with an azimuth bias of −50°).
Frontiers in Earth Science | www.frontiersin.org September 2021 | Volume 9 | Article 716630 Frontiers in Earth Science | www.frontiersin.org September 2021 | Volume 9 | Article 716630 separate from the triggering earthquake, which has been modeled as a deep-crustal earthquake with different rupture processes, and surface wave radiation patterns (Gou et al., 2019;Hua et al., 2019;Kobayashi et al., 2019;Zang et al., 2019). Although we report an event detection, we do not expect significant shaking from radiated energy from the landslide.

Hypothesis II: Local Amplification of Earthquakes from Regional Tectonics
Previous studies have documented the occurrence of earthquakes in the Western African continental region and especially within continental Nigeria (Blundell, 1976;Wright, 1976;Williams and Williams, 1977;Scheidegger and Ajakaiye, 1985;Ajakaiye et al., 1987). These earthquakes have been attributed to 1) the inland continuation of oceanic fracture zones (Blundell, 1976;Wright, 1976) or 2) seismicity associated with faulting along the Cameroon volcanic line (Tabod et al., 1992;Nfomou et al., 2004). The apparent extension of the oceanic fracture zones across the continent-ocean boundary towards the shoreline in the Gulf of Guinea ( Figure 1A; Granot and Dyment, 2015) may explain the clustering of historical seismicity close to the shorelines ( Figure 1A), but may not explain the events reported from areas further inland. We note that historical tremors that were felt in the inland areas (north central part of Nigeria) are distributed across both the failed rifts and areas of the exposed basement on the rift flanks ( Figure 1A). We consider an alternative hypothesis that the series of shakings felt in the study area (Abuja, north central Nigeria) may be related to the local amplification of earthquakes associated with the active tectonics along the offshore oceanic fracture zones (transforms), or along the active Cameroon Volcanic Line in the east. The Abuja area is located within the broad basement complex terrain of Nigeria (i.e., outboard of the Mesozoic failed rifts), where the basement rocks are dominated by Precambrian migmatite gneiss and granitic intrusions (Oyawoye, 1964). The granitic intrusions occur as prominent topographic-highs in the landscape, and the gneissic host rocks occur as topographic-lows where small basins of unconsolidated sediments commonly accumulate. Our estimation of basement undulations in the Abuja area shows that the locations where the tremor was felt are within a small sedimentary basin ( Figure 6G), indicating a vulnerability of the local geology of the Abuja area to seismic amplification. Although the oceanic transforms, coastal areas, and areas near the Cameroon Volcanic Line show historical seismic activity ( Figure 1A), our catalog search did not identify any event along the offshore oceanic transforms or the Cameroon volcanic Line that coincides with the timing of the shaking.

Hypothesis III: Reactivation of Local Basement Fractures and Fault Systems by Post-Rift Crustal Relaxation or Far-Field Tectonic Stresses
Although we report inconclusive evidence for an earthquake in Nigeria, we emphasize that the absence of evidence is not the evidence of absence, considering the poor data quality from the closest station (TORO). As mentioned earlier, the historical earthquakes in Nigeria have been attributed to the inland continuation of oceanic fracture zones through the Precambrian basement shear zones (Blundell, 1976;Wright, 1976;Odeyemi, 1989;Anifowose et al., 2006;Odeyemi, 2006). However, these Precambrian basement shear zones (cyan dotted lines in Figure 1A) show a NNE trend that is markedly distinct from the ENE-to-NE trends of the oceanic fracture zones. Rather, the oceanic fracture zones show better Frontiers in Earth Science | www.frontiersin.org September 2021 | Volume 9 | Article 716630 spatial and azimuthal correlation with the NE-trending Benue Trough (failed Mesozoic rift) than the onshore basement shear zones ( Figure 1A). In addition, there exists no surficial evidence of recent fault scarps or collocated active brittle exploitation along the trend of the basement shear zones. Therefore, we examine the brittle structures in the basement of the Abuja area which is located on the flank of the Benue Trough, with a view of possible local seismic reactivation of pre-existing basement fracture systems. Our analysis of the exposed basement fracture systems ( Figures  6A-C) and deeper structures of the aeromagnetic lineaments ( Figures 6D,E) show common prominent sets of conjugate NE-SW and NW-SE trends. Thus, we interpret these structural trends as the dominant patterns of the basement fracture systems and faults in the Abuja area. We note that the Abuja area is located near the nexus of two Mesozoic rift basins (failed rifts), the NW-trending Bida Basin and NE-trending Benue Trough (Figures 1A, 6C). The basins are composed of basement-rooting structures associated with the opening (Mid Cretaceous) and closing (Late Cretaceous) of the Western and Central African Rift System. We find that the axial trends of these two rift structures are parallel to the most prominent structural trends in the Abuja area ( Figure 6C). This suggests that the pre-existing discontinuities in the Abuja area are most likely inherited from the Mesozoic extensional tectonic deformation. We propose that these basement structures may be reactivated either by far-field tectonic stresses or post-rift relaxation of the crust. This reactivation can be accompanied by seismicity or developed aseismically.
Unfortunately, there exists no published stress field data in the interior of Nigeria (Heidbach et al., 2018). Nevertheless, to provide a first-order assessment, we consider the closest available measurement of the in situ azimuth of 176°for the maximum horizontal compressional stress (SHmax from borehole breakout) at a location in SW Nigeria ( Figure 6C inset; Heidbach et al., 2018). If these available stress data are representative of the current σ1 magnitude and orientation in the Abuja area, the NW-trending fractures are more critically oriented for reactivation (in strike-slip mode) than those in the NE and E-W trends. We recommend future studies to explore the in situ stress field across the region in order to better understand the susceptibility of the inherited structures to seismic reactivation.
Based on the geological considerations presented above, we infer that in the Abuja area, there exist basement discontinuities that could be optimally oriented for seismic reactivation by stress perturbation of the assumed current stress field. On a regional scale of the Western Africa sub-continent, we highlight the occurrence of earthquakes in the Mesozoic rift structures and the Precambrian basement domains (earthquakes in Termit Basin, Benue Trough, Yola-Doba Basin, Ahaggar Massif, and West African Craton in Figure 6A), suggesting that inherited brittle structures in the basement pose important seismic hazards across the region. Overall, we emphasize that the faults and fracture systems of the failed Mesozoic rift basins in the Western Africa region (West African Rift System, Figure 1A) represent critical seismic hazards in the region that may be capable of hosting damaging earthquakes.

Hypothesis IV: Induced Earthquake Triggering by Groundwater Extraction
The widespread subsidence observed in the study area is consistent with those measured in other metro areas around the world and is likely to be associated with fluid (primarily groundwater) extraction (Chaussard et al., 2014;Miller and Shirzaei, 2015;Miller et al., 2017;Miller and Shirzaei, 2019;Herrera-García et al., 2021). A change in near-surface hydrologic loading has the potential to alter the local and regional stress field (Amos et al., 2014;Johnson et al., 2017a;Johnson et al., 2017b;Carlson et al., 2020;Johnson et al., 2020), which can encourage earthquake nucleation and may weakly modulate seismicity (Heki, 2003;Johnson et al., 2017a). In regions where elastic loading maintains a strong periodic signal, the same cyclic pattern is observed in seismic catalogs (Heki, 2003;Ader and Avouac, 2013). Furthermore, stress and pressure in the crust can be altered due to fluid injection and extraction, triggering earthquakes (Wang, 2001;Wang and Kümpel, 2003;Segall, 2010;Ellsworth, 2013;Segall and Lu, 2015;Shirzaei et al., 2016;Keranen and Weingarten, 2018;Kwiatek et al., 2019;Zhai et al., 2019;Zhai et al., 2021). The stress and pressure gain is a function of pumping and injection location, flux, and accumulated volume, while the rate of radial fluid migration is mainly controlled by crustal permeability (Wang, 2001;Segall and Lu, 2015).
Fault failure may occur when shear stress exceeds the fault shear strength for a given effective normal stress. Shear stress can be altered due to non-zero differential stress changes caused by the imparted stress. Also, effective normal stress depends on the magnitude of the stresses and the orientation of the fault in the tectonic stress field. The magnitude of fault-normal stresses can be reduced due to increased pore fluid pressure. However, establishing a threshold for stress and pressure change to trigger an earthquake is not trivial (e.g., Talwani and Acree, 1985). It is often assumed that faults are near critically stressed, if they have not ruptured recently (Townend and Zoback, 2000). Therefore, a small perturbation of the stress field due to the loading effect and fluid diffusion may trigger some earthquakes. Some examples include seasonal modulation of seismicity due to hydrological loading cycles (Christiansen et al., 2007;Johnson et al., 2017a;Carlson et al., 2020), triggering earthquakes due to tides (Wilcock, 2001;Tanaka et al., 2002), and induced seismicity due to pore pressure change by seasonal precipitation or snowmelt (Saar and Manga, 2003;Hainzl et al., 2006;Montgomery-Brown et al., 2019).
Despite scientific evidence from elsewhere, in the Abuja case study, investigating the hypothesis that the events similar to that of September 2018 are of hydrological origin is not straightforward due to the lack of dense hydrological observations (e.g., groundwater levels and stream discharge) and a complete seismic catalog. Therefore, we call for future efforts to generate new observations and develop models that constrain spatiotemporal variations in components of total water storage and establish a link to the local and regional tectonics and seismicity. Such data and models will further enhance the knowledge of water availability and improve the local communitie's resilience to drought in the era of climate change.

Summary: Preferred Hypotheses and Strategies for Further Tests
Our exploration of the limited yet very important geological and geophysical datasets has shed light on the most plausible cause of the shaking experienced in Abuja between September 5 and 7, 2018. Based on the available InSAR analysis, we find the strongest support for the hypotheses that these events were triggered by groundwater extraction (H-IV). Similarly, recently published strain rate modeling from country-wide GNSS network covering the period of 2012-2015 (i.e., before the 2018 tremors; Bawa et al., 2020) shows a localized subsidence zone (negative dilation) that is collocated with a prominent subsidence zone in our analyses (see location b of Figure 7A). This subsidence anomaly is located directly over the Jere Irrigation Project, which was reportedly active at the time of the reported tremors. Following closely, it is the hypothesis that at least one of the felt events could have been triggered by the remote teleseismic event in Japan (H-I). While this idea is plausible, we have no actual data to confirm this hypothesis, because the closest seismic station did not record good enough data to investigate this hypothesis any further (see examples, Han et al., 2017). This again underscores the need for a high-quality geodetic and seismic network in the region. Similarly, a thorough examination of both hypotheses related to local fault reactivation by regional tectonics would require a high quality seismic and geodetic network to evaluate its plausibility for future events (H-II and -III). All the evidence taken together, we point out that, although the Abuja area is located on a stable intraplate region within the continent, transient strain rates, e.g., recharge and discharge of aquifers, can be more important than the background tectonic strain rates when it comes to triggering an earthquake swarm (Sykes, 1978;Wolin et al., 2012;Calais et al., 2016;Gardonio et al., 2018). The Western African crystalline basement is dominated by inherited crustal weaknesses that may then fail due to the transient stress release. Should such events reoccur in the future, further evaluation of our proposed hypotheses would definitely benefit from a renewed geophysical investigation of the Western African region, which is currently grossly understudied and poorly instrumented.

CONCLUSION
We investigated the source of ground shaking reported during September 5-7, 2018, in the Abuja area, central Nigeria. We reviewed previous seismic activity in the region, speculated on how the shaking is related to unique teleseismic events, or may be connected to other alternative explanations, i.e., anthropogenically triggered, or whether regional tectonics and local geology could have made the region more susceptible to triggered fault rupture and amplification of seismic shaking. We explored the spatial and temporal origins of the shaking using seismology and studied the basement structure and surface deformation using aeromagnetic and SAR data. We found the strongest support for seismicity related to anthropogenic groundwater extraction. While other hypotheses cannot be ruled out completely for the case study presented here, we point out that more work is needed to establish a better understanding of the potential connections between inherited basement structure, active regional tectonics, and anthropogenic stress perturbations.

Data and Resources
All the seismic data used in this study are publicly available through USGS and IRIS. The MATLAB code for polarization was originally developed by Matt Haney based on the covariance method and coherency method discussed by Vidale (1986). M_idist (Pawlowicz, 2020) is used to calculate the distance between the station and source location and the corresponding azimuth, and Matlab_TauP based on (Crotwell et al., 1999) is used to calculate the travel time of seismic phases given a list of event parameters (such as origins and magnitude). We use the irisFetch FDSN event web service method (irisFetch.Events) to retrieve seismic data. Information on earthquake location, magnitude, and focal mechanisms are obtained from the NEIC-PDE and USGS catalogs.

DATA AVAILABILITY STATEMENT
Publicly available datasets were analyzed in this study. These data can be found here: https://doi.org/10.17632/W7KKG3G37D.1.

AUTHOR CONTRIBUTIONS
All authors listed have made a substantial, direct, and intellectual contribution to the work and approved it for publication.

ACKNOWLEDGMENTS
The authors acknowledge the use of the BlueHive Linux cluster at the University of Rochester's center for integrated research computing (CIRC). YL acknowledges a summer research grant from URSeismo and funding from the Georgen Institute of Data Science. The authors thank the Nigeria Geological Survey Agency (NGSA) for providing the aeromagnetic data used in this study. They acknowledge many helpful discussions with members of the URSeismo lab: Baowei Liu, Ziqi Zhang, Liam Moser, Trey Brink, and Faner Lin. The facilities of IRIS Data Services, and specifically the IRIS Data Management Center, were used for access to waveforms, related metadata, and/or derived products used in this study. IRIS Data Services are funded through the Seismological Facilities for the Advancement of Geoscience (SAGE) Award of the National Science Foundation under Cooperative Service Agreement EAR-1851048. The satellite ( Figure 6) and aeromagnetic data (Figure 8) used for this research are published on a public repository at (Olugboji, 2020) and included as google KMZ files in the supplementary information files, Datasets S1 and S2. We thank the editor, PA, and two reviewers for their insightful comments which helped to increase the quality of the manuscript.