# Bulletin of the Seismological Society of America

## Abstract

To date little is known about seismic hazard in Eritrea, despite its location in a volcanically and tectonically active region, and the gathering pace of major infrastructure projects. In response, we report the findings of a comprehensive probabilistic seismic‐hazard assessment for Eritrea and adjacent areas. Seismic source and ground‐motion models are constructed separately; we use an adaptive spatiotemporal smoothing method to map expected patterns of seismicity. To construct a consistent earthquake catalog from different data sets, we use orthogonal regression to convert and unify different magnitude scales. A sensitivity analysis of the different input parameters helps constrain them and disaggregation of site‐specific hazard estimates yields insights into the relative contribution from seismic sources of different magnitudes and distances. The results highlight seismic hazard in proximity to the Red Sea, Gulf of Aden, Afar depression, and along the boundaries of the Danakil microplate. We estimate a 10% chance over 50 years of observing pseudospectral accelerations (PSAs) at 0.2 s exceeding 0.16*g* in the port city of Massawa (population ∼32,000) and the town of Bada (population ∼4000). For the capital, Asmara (population ∼520,000), we calculate a PSA of 0.11*g* at 0.2 s. Compared with previous studies, our results provide greater spatial resolution, use more recent ground‐motion models, and benefit from a smoothed seismicity method. Our aims are to stimulate further studies and contribute to the safe development of the region in light of its exposure to seismic hazards.

## Introduction

Eritrea is located at the northernmost extent of the East African rift system (EARS). As such, it lies within a tectonically active region bounded by the Red Sea to the northeast and the Danakil depression to the southwest. It is prone to earthquakes and volcanic eruptions synonymous with the rifting process; historical records attest to large earthquakes up to surface‐wave magnitude (*M*_{s}) 6.8 (Gouin, 1979) occurring within the country’s borders. Seismic‐hazard assessments of this region have been previously conducted (Gouin, 1976; Kebede and Asfaw, 1996; Midzi *et al.*, 1999; Hagos *et al.*, 2006). Gouin (1976) produced peak ground acceleration (PGA) maps for Ethiopia (which at that time included Eritrea). This was followed by a seismic‐hazard assessment for Ethiopia and neighboring countries (Kebede and Asfaw, 1996). As part of the Global Seismic Hazard Assessment Project (GSHAP), a standard probabilistic seismic‐hazard assessment (PSHA) was conducted over a much wider area in eastern and southern Africa, extending from latitudes 40° S to 20° N (Midzi *et al.*, 1999). These studies lacked detail at the scale of a country such as Eritrea (e.g., Midzi *et al.*, 1999) and used older and less accurate ground‐motion prediction equations (GMPEs) that tend to overestimate seismic hazard (e.g., Hagos *et al.*, 2006). Furthermore, recent deployments of broadband seismic networks have improved our knowledge of the seismicity of the region (Goitom *et al.*, 2015).

Based on the earthquake magnitudes recorded to date, Eritrea and the surrounding area can be considered as a moderate seismicity region. PSHA studies in such regions are challenging owing to the paucity and incompleteness of data, particularly when local stations are lacking. However, seismic‐hazard assessment has become increasingly important in the context of sustainable development (Kracke and Heinrich, 2004; Hamouda, 2011; Goda *et al.*, 2013; Lam *et al.*, 2015).

In this article, we assess the seismic hazard of Eritrea by compiling an earthquake catalog from 1400 to present and using an adaptive spatiotemporal smoothing approach to produce annual seismicity rates. We show that several major towns are exposed to significant seismic hazard, highlighting the need to consider earthquake hazards in construction projects. Other innovations of our study include (1) the use of a general orthogonal regression to convert and unify different scales of magnitudes available to us from different data sources; and (2) more recent applicable GMPEs from the literature (there are no GMPEs specific to the region). Finally, we perform a sensitivity analysis to examine the effects of different input parameters and to understand the potential errors in our assessment.

## Tectonic Setting

The EARS is a classic example of a continental rift. The tectonics of the study area are governed by the divergence of the Arabian, Nubian, and Somali plates (Fig. 1)—which meet at the Afar triple junction, close to the international borders of Eritrea, Ethiopia, and Djibouti.

The Red Sea rift was formed by the separation of the Arabian and Nubian plates in the Early Miocene (Steckler *et al.*, 1988; McClusky *et al.*, 2003). The EARS marks the boundary between the Nubian plate to the west and the Somali plate to the east; both plates are part of what was formerly known as the African plate (Chu and Gordon, 1998). The separation of the Arabian plate from the Somali plate formed the Gulf of Aden. Kinematic models constrained by Global Positioning System (GPS) data suggest ∼15 mm/yr of N35°E directed opening across the western Gulf of Aden rift (e.g., Jestin *et al.*, 1994). Current opening across the kinematically complex southern Red Sea rift is well constrained with a relatively high density of GPS measurements (e.g., ArRajehi *et al.*, 2010; McClusky *et al.*, 2010). These data show that south of ∼17° N, the rift bifurcates into two branches: the main Red Sea and the subaerial Red Sea rift in Afar (Danakil depression; Fig. 1). Partitioning of extension between rift branches varies along strike. North of ∼16°, all the extension is accommodated in the main Red Sea rift, spreading at ∼15 mm/yr. Moving south of 16° N, the extension is progressively accommodated in the Afar depression reaching ∼20 mm/yr at ∼13° N (McClusky *et al.*, 2010).

Although we focus on the hazards in Eritrea, our study area is located between latitudes 10° N and 20° N and longitudes 36° E and 44° E, which encompasses Eritrea, parts of Ethiopia, Sudan, and Djibouti (Fig. 1). In Eritrea, four main physiographic units associated with various tectonic structures can be identified: the plateau, the rift margin, the Afar rift, and the Red Sea rift (Fig. 1; Ogubazghi and Woldehaimanot, 1998). The plateau covers the central and western parts of Eritrea and consists largely of structures of late Precambrian age (Drury *et al.*, 1994). The rift margin is defined as the area between the plateau and the Red Sea margin and the Afar depression. It is characterized by graben and salt planes in the north and south, respectively. The Afar depression is located in southeastern Eritrea and northeastern Ethiopia. It is surrounded by the plateau, the Red Sea, the Somali plate, and the Danakil horst. Chu and Gordon (1998) showed that the seismicity along the Red Sea axis branches into two directions, southeast and south‐southwest at 17.5° N. A similar pattern is also apparent in our seismicity map (Fig. 2). The aseismic block surrounded by the Red Sea and a line of epicenters trending west of Gulf of Tadjurah is the Danakil microplate (Chu and Gordon, 1998).

## Historical Seismicity

The earliest recorded earthquake in Eritrea occurred in 1400 near the Dubbi volcano in southeast Eritrea (Gouin, 1979). It was accompanied by a volcanic explosion that was observed and reported by sailors on the Red Sea and lasted for more than one week (Gouin, 1979). Since then, a number of significant earthquakes have been reported in this region, mostly linked with volcanic activity within the Bidu volcanic massif, a northeast–southwest‐trending chain of volcanoes (Wiart and Oppenheimer, 2005). In May 1861, an *M*_{s} 5.5 earthquake was recorded and associated with the Dubbi volcano eruption (Gouin, 1979). This stands as the largest documented historical eruption on the African continent (Wiart and Oppenheimer, 2000).

More recently, in March 2011, an earthquake of *M*_{L} 4.8 was recorded around Siroru, a village within the Nabro volcanic caldera, roughly 25 km south of Dubbi volcano. This earthquake caused significant building damage and landslides, killing a number of goats (Ogubazghi *et al.*, 2014). This was followed by a large eruption of the Nabro volcano in June 2011 (Goitom *et al.*, 2015).

Not all earthquake occurrences in Eritrea are related to volcanic eruptions. Along the rift axis in the northern Red Sea region, a series of earthquakes were reported around Alid volcano in 1901 and 1902 (Gouin, 1979) and also around Bada in 1993 (Ogubazghi *et al.*, 2004; Ogubazghi and Goitom, 2015). Although these were unrelated to eruptive activity, they were interpreted as being due to magmatic processes in the crust. Massawa, the main port city in Eritrea, is situated on large border faults that demarcate the edge of the Red Sea rift. In 1884, it was severely damaged by an *M*_{s} 5.9 earthquake and has been affected by a number of subsequent moderate earthquakes in 1912, 1913, 1921 (Gouin, 1979), and 2002 (Ogubazghi and Goitom, 2015). Earthquakes were also felt in 1875 and 1876 around the city of Keren, a large town close to the border faults north of Massawa.

As well as earthquakes within Eritrea presenting a potential hazard, earthquakes occur outside its borders that could induce hazardous ground shaking within Eritrea. For example, the Afar depression to the south of Eritrea has a well‐characterized history of seismic and volcanic activity, which has been studied through a number of collaborative projects (e.g., the Ethiopia Afar Geoscientific Lithospheric Experiment and Afar Consortium projects: Ayele, Jacques, *et al.*, 2007; Ayele, Stuart, *et al.*, 2007; Ebinger *et al.*, 2008; Belachew *et al.*, 2011; Hammond *et al.*, 2011; Keir *et al.*, 2013). Major seismic events include an *M*_{L} 6.3 earthquake in Serdo (Ethiopia) in 1969, an *M*_{L} 6.2 earthquake in the Dobi graben (Ethiopia) in 1989 (Ogubazghi and Woldehaimanot, 1998), and a series of dike injections in 2005–2010 that induced seismicity in Mandar‐Harraro (Ethiopia) (e.g., Ayele *et al.*, 2009).

## Methodology

Our methodology follows the framework of PSHA (Cornell, 1968), with particular model choices informed by recent developments in seismic source and ground‐motion modeling and constrained by the availability and quality of data in our study region.

PSHA requires a seismic source model and a ground‐motion model (e.g., Ambraseys *et al.*, 1996; Chiou and Youngs, 2008; Akkar and Bommer, 2010). Depending on the availability and quality of earthquake data, the seismic source model is often specified with estimated seismicity rates in distinct tectonic and seismological zones delineated by experts. These zones may be augmented with more specific earthquake source areas, often associated with mapped faults, which are assumed to rupture with characteristic magnitudes (Hodge *et al.*, 2015).

Smoothed seismicity models have become a popular alternative to tectonic zonation for estimating the spatial distribution of seismic rates, see for example, the U.S. National Seismic Hazard Map (Petersen *et al.*, 2014) and the European hazard map (Woessner *et al.*, 2015). Smoothed seismicity models have shown promising predictive skill in prospective, out‐of‐sample testing (Schorlemmer *et al.*, 2010; Werner *et al.*, 2010; Zechar *et al.*, 2013). They are also more objective and reproducible than models using tectonic zones. In this study, we use an adaptive smoothing method (e.g., Stock and Smith, 2002; Helmstetter *et al.*, 2007; Helmstetter and Werner, 2012, 2014), which is increasingly used to calculate seismicity rates for official government agency hazard models (e.g., Felzer, 2013; Field *et al.*, 2014).

### Catalog Compilation

The starting point for PSHA is a complete earthquake catalog that includes preinstrumental, early instrumental, and modern instrumental data. The study area is characterized by shallow to medium depth and moderate size earthquakes governed by the tectonic activities of the Red Sea, Gulf of Aden, and East African rift (Ayele and Kulhanek, 1997; Hagos *et al.*, 2006).

As a first step to the seismic‐hazard assessment of the region, all available catalogs were assembled. We sourced the following to compile the seismic catalog of the region.

Gouin (1979) was used for earthquakes covering the time span from 1400 to 1964. It reports historical earthquakes of the Horn of Africa compiled by the author through consulting libraries, churches, newspapers, and different scholars of the region. It lists and describes earthquakes in terms of location, timing, magnitude/intensity, and felt localities. This remarkable work is the most comprehensive compilation of historical earthquakes for the region. In some cases, however, care needs to be taken as epicentral locations and depths are not accurate and most of the events pre‐1900 lack magnitudes or intensities.

We use the reviewed version of the International Seismological Centre (ISC) earthquake catalog that covers the time span from 1964 to 2012 (ISC, 2012).

We use a local catalog that covers the period from 1999 to 2004 and that was compiled by the University of Asmara and the Eritrea Institute of Technology through locating events from short‐period seismic stations deployed in the country.

An additional data set was acquired from a regional organization, the Eastern and Southern African Regional Seismological Working Group (ESARSWG), and covers the period 1993 to 2012. This working group includes nine countries of eastern and southern Africa and has prepared an earthquake catalog of the region through data sharing at annual workshops.

We also included the earthquakes that we located from six broadband stations deployed in Eritrea together with 10 stations in Ethiopia spanning the period from June 2011 to October 2012 (Illsley‐Kemp

*et al.*, 2017).

We adopt the criteria set by Nasir *et al.* (2013) to remove duplicate events in which events that are located within 10 km and within occurrence time of 2 min are considered as duplicates. Whenever an event is listed in more than one catalog an order of priority—local, regional, ISC, and Gouin—was considered for retaining the event in the catalog. Based on this, 648 events were found to be duplicates and removed from the catalog. The resulting catalog contains 5896 events between 1400 and 2012 (Fig. 2a); and Figure 2b shows the temporal variation of magnitudes for 1900–2012 for all catalogs.

To address any bias in the composite catalogs by including regional and local earthquake records, we use two different catalogs to assess the hazard. Model 1 uses the catalog compiled from the ISC database that only covers the duration between 1964 and 2012 and model 2 includes the local catalogs from Eritrea and ESARSWG for the periods 1999–2004 and 1993–2012, respectively. A third model, which includes the historical earthquakes in addition to model 2, is also evaluated.

### Magnitude Conversion Using Orthogonal Regression

The magnitudes in the composite catalog that have been compiled from different catalogs are of different scales. Notably there are two magnitude scales: local (*M*_{L}) and body magnitudes (*m*_{b}). Although the ISC catalog reports both local and body magnitudes, the local and the ESARSWG catalogs only provide local magnitudes. In PSHA, moment magnitude (*M*_{w}) is commonly preferred because it does not saturate for higher magnitudes and is also related to a physical property of the source.

There are two methods described in the literature that are used to convert *m*_{b} to *M*_{w}: linear regression and orthogonal regression. The PSHA applied in this region so far used a linear regression method to convert between magnitude scales (Kebede and Asfaw, 1996; Midzi *et al.*, 1999; Hagos *et al.*, 2006). Most studies have used the global relation developed by Scordilis (2006). However, some authors (e.g., Castellaro *et al.*, 2006; Das *et al.*, 2013) have shown that linear regression has limitations. For example, it assumes that the error of the independent variable is negligible. Instead, Das *et al.* (2013) applied an orthogonal regression, which takes into account the effects of measurement error in both variables, to convert *m*_{b} to *M*_{w} in four regions (Japan, Mexico, Indian Himalaya, and Peninsular India) and showed that the orthogonal regression method resulted in lower average differences and lower standard deviations.

To implement a similar orthogonal regression approach for our study area, 278 events have been extracted from the ISC (for *m*_{b}); and from the Global Earthquake Model (Storchak *et al.*, 2013) and the Global Centroid Moment Tensor (Dziewonsky *et al.*, 1981; Ekstrom *et al.*, 2012) catalogs (for *M*_{w}), and used in both linear and orthogonal regressions. A total of 250 events are used to develop the relationships, and the remaining 28 events are used to check the validity of the derived relationships. The procedure for orthogonal regression is described in detail in the literature (e.g., Castellaro *et al.*, 2006; Das *et al.*, 2011, 2012). The two magnitudes *M*_{w} and *m*_{b} are first fitted in an orthogonal manner, yielding the relation: (1)We use *m*_{b(proxy)} instead of *m*_{b} in equation (1) because of the criteria of minimization of orthogonal residuals to a best‐fit line. We need to obtain *m*_{b(proxy)} to predict *M*_{w} (Das *et al.*, 2012, 2013). *m*_{b(proxy)} is obtained from the abscissa by projecting the points (*m*_{b}, *M*_{w}) orthogonally to the fitted line. Then, a linear regression is applied between the new points *m*_{b(proxy)} and *m*_{b} to find a relation between the two: (2)Inserting equation (2) into (1), we find a relation between *M*_{w} and *m*_{b},(3)which is used to convert *m*_{b} to *M*_{w}. A similar approach has been used to develop a relation between *M*_{L} and *M*_{w}, (4)We compare the results of the orthogonal regression with a linear regression, (5)and the global relation by Scordilis (2006), (6)Our linear regression (5) is similar to the global relation (6). However, the orthogonal relation (3) is different from these two. We compare the fit using the mean absolute difference (MAD) between *M*_{w} and *m*_{b} from the 28 events. The MADs for orthogonal, linear, and Scordilis relations are found to be 0.123, 0.128, and 0.127, respectively. The results of the three methods for the 28 events are shown in Figure 3. Comparing the three methods in terms of magnitude ranges, we observe that the differences are smaller at higher magnitudes (*M*_{w}>5.0) and more pronounced for *M*_{w}<5.0. Although the differences are not large, the orthogonal regression results in the smallest MAD; this is in agreement with results from studies in other regions (e.g., Das *et al.*, 2013). Based on these results, we use orthogonal regression to unify and convert the body and local magnitudes to moment magnitude using equations (3) and (4), respectively.

A number of completeness analysis methods are available (e.g., Mignan and Woessner, 2012; Nasir *et al.*, 2013). We used a maximum curvature approach (Wiemer and Wyss, 2000) to determine the magnitude completeness of our catalog (Fig. 4). The histogram reveals the maximum number of earthquakes for *M*_{w} 4.0 and hence it may be considered that the catalog is complete for *M*_{w}≥4.0. The bimodal features can be seen in the histogram because of inclusion of more data of smaller magnitudes from local stations after 1993. Conservatively, we considered that our catalog is complete above magnitude 4.5 since 1964.

### Recurrence Parameters (*b*‐Values)

The *b*‐value that describes the relationship between the numbers of small to large earthquakes is one of the input parameters in the study of seismic hazards. The well‐known Gutenberg–Richter (GR) relationship is (7)in which *N* is the number of events with magnitudes equal to or greater than *m*, and *a* and *b* are constants. It describes the recurrence relation and was used to compute *b*‐values for the study area.

We used a maximum‐likelihood approach (Williams and Calvez, 2012) to compute a uniform *b*‐value for the whole region. In this process, the minimum magnitude completeness and error limits (upper and lower) are calculated from the input magnitudes.

The *b*‐values computed are 1.09 (with a 95% confidence interval from 1.00 to 1.17) and 1.10 (from 0.99 to 1.18) for models 1 and 2, respectively (Fig. 5). The *b*‐values for both models are similar (1.1) and found to be within the range found in earlier studies of the region (Kebede and Kulhanek, 1994; Ayele and Kulhanek, 1997; Jacques *et al.*, 1999; Hofstetter and Beyth, 2003; Hamlyn *et al.*, 2014; Illsley‐Kemp *et al.*, 2017). Figure 5 also shows the annual seismicity rates between 1964 and 2012 for both models. The annual earthquake rates of *M*_{w}≥4.5 are ∼18 and 19 for models 1 and 2, respectively.

### Adaptive Smoothing

The PSHA studies conducted so far in East Africa are based on seismic zoning and spatially smoothed seismicity to produce seismicity rates (Midzi *et al.*, 1999; Hagos *et al.*, 2006). Hodge *et al.* (2015) illustrated the use of geomorphology and geodesy together with seismicity to estimate hazards in the Malawi rift. We use an adaptive spatiotemporal smoothing method to estimate seismicity rate for the region. This method was developed by Helmstetter and Werner (2012) and applied in California. It avoids the need to divide the region into source zones. It also accounts for the spatiotemporal variability of the catalog completeness by correcting for potentially missing earthquakes assuming a GR relationship.

The catalog is divided into two parts: learning and target catalogs. The events in the learning catalog are used to construct a predictive earthquake density. They form the centers of kernels with adaptive widths (in space and time). Spatial adaptive bandwidths are proportional to the distance to the *k*th nearest neighbor, in which *k* is a parameter. This can be generalized to space–time (for details, see Helmstetter and Werner, 2012, and references therein). The earthquakes in the learning catalog are used to construct a spatiotemporal probability density as the normalized sum of the kernels. Parameters (such as *k*) are optimized by reducing the misfit between the constructed earthquake density (the model) and the events in the target catalog. The misfit, or inversely the predictive skill, of the density is calculated as the joint log likelihood of the target earthquakes given the constructed density. Parameters are estimated by maximizing the joint log likelihood.

The density is calculated in 0.1° by 0.1° cells every *D*_{t} days using the space–time kernels. The contribution to this density is greatest from nearby and recent earthquakes, and the contribution decreases for more distant and less recent events according to the smoothness of the kernels. The estimated seismicity rate in each cell effectively becomes a time series (rate every *D*_{t} days) that tracks clustering. To estimate long‐term seismicity rates, temporal clustering effects are removed by taking the median of the estimated seismicity rate time series in each cell. This amounts to a nonparametric method for removing foreshocks and aftershocks that is more objective and less sensitive to arbitrary parameters than standard deterministic declustering algorithms (Helmstetter and Werner, 2012).

Based on this method, we compute multiple models by varying the input parameters and data sets. We optimize the parameters for each model and analyze the likelihood scores of the models with respect to a benchmark (reference) model to infer optimal model choices. The reference model is a spatially uniform density. We compute the probability gain, defined as the probability gain per target earthquake relative to a spatially uniform Poisson model, of each model over the benchmark as a measure of predictive skill. The probability gain is given by the following equation: (8)(equation 12 in Helmstetter and Werner, 2012), in which *N*_{u}=*N*_{t}/*N*_{c}, *N*_{t} is the total number of targets, *N*_{c} is the total number of cells in the testing area, and *N*_{p} is the predicted rate in the bin.

### Monte Carlo Simulation of Synthetic Catalogs

We used a Monte Carlo approach to produce synthetic catalogs. Because the duration of simulation has significant influence on the PSHA results, preliminary investigations were conducted by varying the simulation duration from 100,000 to 5,000,000 years. It was found that results based on synthetic catalogs for 1 million years produce stable seismic‐hazard estimates, whereas increase in the simulation duration leads to significantly high computational effort, especially when seismic‐hazard contour maps are developed. Based on these preliminary investigations, we decided to present the PSHA results based on the synthetic catalog having a 1‐million‐year duration. Simulations were based on random selection of epicentral locations from the spatial probability density, five depths (10, 15, 20, 25, and 30 km with equal probability), and magnitudes from the GR relationship (with maximum magnitudes of 7, 7.25, and 7.5 with weights 0.25, 0.5, and 0.25, respectively, and *b*‐values of 0.9, 1.0, and 1.1 with weights 0.25, 0.5, and 0.25). In our seismic‐hazard calculations, a hypothetical fault geometry is generated for each event to account for the effects of finite‐fault dimensions on the ground‐motion prediction and is commonly determined through fault‐rupture parameters (Atkinson and Goda, 2011). Because detailed fault studies of the area are not yet available, we assumed a dip angle of 90° by taking the source parameters that have been determined based on empirical scaling relationships by Wells and Coppersmith (1994). We acknowledge that these scaling relationships were developed more than two decades ago and are obsolete, and newer relationships have been developed (e.g., Leonard, 2010; Yen and Ma, 2011). However, we used them for comparison purposes and to investigate the use of probabilistic conversion models of different distance measures in PSHA (Scherbaum *et al.*, 2004; Goda *et al.*, 2010) to reduce computational efforts.

### Ground‐Motion Model

GMPEs are important model components for PSHA. A GMPE characterizes the propagation of seismic waves as a function of magnitude and source‐to‐site distance. Strong ground motion data are needed to develop a GMPE for a region. A suitable GMPE for the region has not yet been developed for use in PSHA (Midzi *et al.*, 1999; Hagos *et al.*, 2006; Hodge *et al.*, 2015). The option therefore is to use available GMPEs that are developed for other regions with a similar tectonic environment. Hagos *et al.* (2006) used the GMPE developed by Ambraseys *et al.* (1996), but this relation tends to overestimate hazards (Ambraseys *et al.*, 2005) and more recent GMPEs have been developed to address this. We adopted four GMPEs: Campbell and Bozorgnia (2008), Chiou and Youngs (2008), Boore and Atkinson (2008), and Akkar and Bommer (2010) which were developed for active tectonic regions characterized by shallow crustal faulting in western United States, Europe, Mediterranean regions, and the Middle East. Although these are the models used in many PSHA studies, they are not the newest ones; with the availability of additional information in the future, the results can be improved using more appropriate GMPEs. We tested different weights for these GMPEs, the results of which are detailed in the Sensitivity Analysis section.

### Probabilistic Seismic‐Hazard Analysis

PSHA is one of the widely used methods for seismic‐hazard assessment. It was developed by Cornell (1968) and later extended by various authors and applied to different seismically active regions. It integrates spatiotemporal seismicity and earthquake occurrences, earthquake size scaling and GMPEs to evaluate the expected ground motions. By considering all possible scenarios, it quantifies the annual probability of exceeding ground‐motion levels at a given location. The approach produces contour maps of the expected ground shaking over a specified time period in a region. Ground shaking is often expressed by PGA or another quantity of interest. PSHA can also produce pseudospectral acceleration (PSA) curves at specific sites within the area, which are additional indicators of hazard levels.

The generalized equation of probabilistic hazard is given by (9)in which *v*(*z*) is the mean annual rate of exceedance of a ground‐motion level *z* at a site of interest, *λ*_{i} is the mean frequency of occurrence of earthquakes exceeding *M*_{min} at a source *i*, *f*_{i}(*m*) and *f*_{i}(*r*) are probability density functions of magnitude and distance, respectively, and *P*(*Z*>*z*|*m*,*r*) is given by a GMPE. The summation is over all *n* sources.

The code we developed to compute the seismic hazard consists of generating synthetic catalogs based on the smoothed seismicity rate and produces PGA and PSA contour maps as well as seismic‐hazard curves and uniform‐hazard spectra at different return periods (e.g., 500, 2500, and 5000 years) for selected sites. We also disaggregated the hazard for one of the cities (Massawa) at a return period of 500 years. In addition to the *b*‐values and the GMPE, the input parameters to the PSHA include the completeness magnitude above which the catalog is complete (*M*_{c}), the minimum magnitude below which no structural damage is expected (*M*_{min}), and the maximum magnitude—which is the maximum expected magnitude in the region (*M*_{max}). These parameters were set to 4.5, 4.0, and 7.5, respectively, by evaluating the catalogs and the vulnerability of buildings in Eritrea. Knowledge of local seismic site conditions measured through *V*_{S30} is also an important factor in PSHA. This is a reference velocity for the upper 30 m and is required in the calculation of the seismic hazard. We took a value of 760 m/s as a reference velocity. We have not modeled active faults because they have not yet been comprehensively mapped in Eritrea. As such data become available, they should be incorporated into the PSHA.

## Results and Discussion

In the production of the smoothed seismicity rates, a number of models with various input parameters were tested: time step (*D*_{t}), minimum magnitude for generation of model (*M*_{d}), minimum magnitude for target events (*M*_{t}), and duration of target catalog (*T*_{t}). Each parameter was tested with different values while keeping others fixed. Table 1 and Figure 6 show the probability gain for the different parameters. Varying the time step *D*_{t} (10, 20, 100, and 500 days) yields similar probability gains. When the duration for the target catalog is short, the probability gain is high indicating that shorter duration of target catalog (this means that the duration of the learning catalog is increased) gives rise to better results. There are no unequivocal trends in the variation of the gain with magnitude thresholds. However, in general the gain is high when the minimum magnitude of the learning catalog is lower and the magnitude for the target catalog is higher.

### Comparison of Models

Before exploring the models fully, we compare the three models that are based on different input catalogs (model 1: ISC [1964–2012]; model 2: local and ISC catalogs [1964–2012]; and model 3: historical [1900–1964], local, and ISC catalogs [1964–2012]). We contrast their predicted hazard levels in Massawa (Fig. 7). The PSAs expected to be exceeded with 10% probability over 50 years are similar for models 1 and 2 across all considered vibration periods. On the other hand, model 3, which also uses historical earthquakes as model input, estimates significantly higher accelerations than models 1 and 2 at periods below 1 s. This is because some large historical earthquakes occurred around Massawa. This fact has to be kept in mind when interpreting the results of model 1 or 2.

### Smoothed Seismicity Rates

The seismicity rates for models 1, 2, and 3 are shown in Figure 8. Models 1 and 2 reveal a similar pattern with minor variations, due to the small number of earthquakes included from the local catalogs in model 2. However, model 3, which includes historical earthquakes from 1900, shows a different pattern. We have less confidence in model 3 because of the historical catalog’s incompleteness. The number of earthquakes in the longer time span of the historical catalog is small (207 earthquakes in 64 years); many low‐magnitude events must have been missed—either there is no information or they were not felt by people (e.g., earthquakes that occurred in the sea). The historical catalog has a higher completeness magnitude and contains few Red Sea earthquakes. As a result, the median values over a longer time span of the catalog remain low along the Red Sea rift resulting in a different pattern for model 3. The temporal variations in the spatial completeness magnitude pose a challenge for seismicity rate estimations, which is beyond the scope of our present work. However, as shown in Figure 8, we believe that excluding the historical earthquakes may underestimate the hazard levels in some sites and care should be taken when interpreting the seismicity rates of model 1 or 2.

All of the western part of the plateau, which is roughly west of the rift margin, is characterized by low seismicity rates. In contrast, the Red Sea rift, the Gulf of Aden, and the northern tip of East African rift show high seismicity rates with the areas around Massawa, Zula, and Alid‐Bada in Eritrea being characterized by particularly elevated seismicity rates. There is a very good correlation between the high seismicity rate and the location of the Danakil microplate (McClusky *et al.*, 2003). Nabro volcano, where seismic activity was high in 2011 (Goitom *et al.*, 2015), is not clearly marked by a high expected seismicity rate; this reflects the rate estimation, which uses the median of the seismic rate time series estimated every *D*_{t} days over the learning catalog. Despite the recent seismicity, the median rate remains low and our results have underestimated the hazard due to the volcanic eruption.

Our results were compared with those of Hagos *et al.* (2006). Although the general trends look similar in that both models (their fig. 4a and our Fig. 8a) show highest seismicity rates around the Gulf of Aden and Red Sea rift, differences are evident. The Red Sea rift is more detailed in our model and the Massawa area is characterized by higher seismicity rates compared with their model. We believe the differences arise because (1) we use adaptive kernels in space and time, whereas they use fixed‐width kernels in space only, (2) we use a longer time span and more events for the input catalog, compared with theirs, and (3) we use the median value in each cell, whereas they count the number of earthquakes in each cell.

### Seismic‐Hazard Maps

Following the procedures outlined in the Methodology section, we calculate the hazard according to models 1 and 2 to produce hazard maps (Fig. 9). The hazard levels generally follow the spatial trends of the seismicity rates. At the 500‐year return period (Fig. 9a, which is equivalent to 10% exceedance level in 50 years), hazard levels greater than 0.1*g* are observed along the plate boundaries: along the Red Sea rift, the Gulf of Aden, and the Afar depression. A band of slightly elevated hazard level of greater than 0.07*g* extends toward the west and northwest of the Gulf of Aden and bends toward the northeast from Massawa to the Red Sea rift. This hazard level roughly follows the western and northern boundaries of the Danakil microplate along the Alid‐Bada line and through Massawa toward the Red Sea rift. The cities of Massawa and Bada fall within this band and Asmara (the capital city of Eritrea), Tio, and Asseb (a port city) are on the margins of this band. The hazard levels increase with increasing return periods but the spatial pattern remains similar.

### Site‐Specific Seismic‐Hazard Curves, Uniform‐Hazard Spectra, and Disaggregation

We selected eight sites: Asmara, Massawa, Asseb, Bada, Tio, Afambo, Nakfa, and Keren (population shown in Table 2) within Eritrea and estimated site‐specific hazard levels for model 1 (the sites are indicated in Figs. 1 and 9). In Figure 10, we show resulting PSA curves against vibration periods for return periods of 500, 2500, and 5000 years. The hazard levels of all the sites increase with higher return periods. Taking Massawa as an example, the PSA at 0.2 s is 0.16*g*, 0.37*g*, and 0.50*g* for return periods of 500, 2500, and 5000 years, respectively. For all return periods the highest PSA is observed in Massawa and Bada, and the lowest in Nakfa. Although the population of Bada is small, as seen from Table 2, with the planned mining development around the area (see Data and Resources), the population is expected to increase soon. With the possible start of potash mining activity around Bada and its proximity to the rift, the effect of induced seismicity should be considered and monitoring policy should be in place during the mining activity.

As found in similar studies, the highest PSA is observed at the vibration period of 0.2 s and at longer vibration periods the PSA is decreasing. At the 500‐year return period, the PSA at 0.2 s are 0.11*g*, 0.16*g*, 0.16*g*, 0.13*g*, 0.11*g*, 0.11*g*, 0.06*g*, and 0.08*g* for Asmara, Massawa, Bada, Tio, Asseb, Afambo, Nakfa, and Keren, respectively. We produced seismic‐hazard curves as a function of peak spectral acceleration for all eight cities but only a comparison of the hazard levels at a vibration period of 0.2 s for Massawa, Tio, Asmara, and Keren is presented in Figure 11 for the sake of brevity. Among these, Massawa has the highest hazard level followed by Tio, Asmara, and Keren. The hazard level decreases when the distance of the site from the rift increases, indicating the rift is the main source of seismic hazard in Eritrea.

To determine the greatest contributing seismic sources to the hazard in Massawa, we disaggregate the 500‐year return period hazard in terms of distance and magnitude at vibration periods of 0.1, 0.2, 1, and 2 s. The results are shown in Figure 12. At vibration periods of 0.1 and 0.2 s in which the PSA is the highest, earthquakes of magnitudes ranging between 4.5 and 6.0 at distances less than 100 km are the major contributors to the hazard levels in Massawa. For example, the earthquakes with magnitude around 7 at distances greater than 150 km (at 0.1 s) contribute little to the hazard level in Massawa, as indicated by the small amplitude in the probability mass function (Fig. 12a). The scale bar in Figure 12 shows , which is a measure of deviation from the median of the predicted PSA for a given pair of magnitude and distance. With increasing vibration periods, the contributions of large earthquakes at larger distances become higher. At vibration periods 1 and 2 s, the contributing earthquakes are distributed more broadly within a distance of 300 km (Fig. 12c,d). These findings have direct implications for building design; when multistory buildings with long natural vibration periods are to be constructed in Massawa, potential earthquake sources up to 300 km away should be considered.

### Sensitivity Analysis

In PSHA, there are a number of input parameters with various degrees of uncertainties and in most cases a logic‐tree approach is used (Kulkarni *et al.*, 1984; Stucchi *et al.*, 2011; Marzocchi *et al.*, 2015). To assess the effects of the choice of the different parameters on the hazard calculation, we perform a sensitivity analysis of the hazard levels at the return period of 500 and 5000 years for Massawa. The parameters selected for the sensitivity analysis were maximum magnitude, minimum magnitude, maximum distance, GMPE, and input catalogs.

Figures 13 and 14 show the results of this analysis for the four parameters for return periods of 500 and 5000 years, respectively. The choice of maximum magnitude (6.5, 7.0, 7.5, and 8.0) does not strongly influence the hazard level of Massawa (Figs. 13a and 14a). The hazard levels are identical for maximum magnitudes 7.5 and 8.0 because their occurrence is so unlikely that they do not contribute to the shaking over the considered time periods. At a vibration period of 0.2 s, the hazard levels increase only very slightly, by 0.5% and 0.06%, when the maximum magnitude is changed from 6.5 to 7.0 and to 7.5, respectively. In our preferred hazard map, we use *M*_{w} 7.5 as the maximum magnitude, which is slightly higher than the observed maximum magnitude of the catalog (*M*_{w} 7.1). This value is consistent with a global study by Bird and Kagan (2004), who estimated a corner magnitude of 7.64 for continental rift boundaries. A corner (rather than maximum) magnitude defines a crossover point in the upper tail of the magnitude distribution from the exponential GR law to a much faster decay.

When the minimum magnitude changes from 4.5 to 5.0, the hazard decreases from 0.16*g* to 0.13*g* at the vibration periods of 0.2 s (i.e., by about 15%; Fig. 13b). Therefore, taking a higher minimum magnitude tends to underestimate the hazard level in Massawa at 500‐ and 5000‐year return periods, which has also been reported in similar hazard studies. Hence, we recommend and use a minimum magnitude of 4.0, which is slightly less than the completeness magnitude, that is, we extrapolate the observed seismicity rates to lower magnitudes using GR scaling.

There is very little influence of the choice of the maximum distance on the results. When, at the 500‐year return period, we change the maximum distance from 100 to 200 km the hazard level in Massawa increases by 4%. When the maximum distance increases from 200 to 300 km, the hazard increases by 4.5%, and it increases by 2.8% when the distance is changed from 300 to 400 km. This shows that the change in the hazard level due to the change in the maximum distance at both return periods (500 and 5000 years) is negligible. Hence, 300 km is used as the maximum distance for the calculation of hazard in the region and this is consistent with the disaggregation results discussed above.

We tested four GMPEs by employing each individually and by combining all with equal weights. The results are shown in Figures 13d and 14d. At a vibration period of 0.2 s, the Chiou and Youngs (2008) model yielded the highest hazard level of 0.180*g* followed by the models by Boore and Atkinson (2008), Campbell and Bozorgnia (2008), and Akkar and Bommer (2010). Table 3 lists PSAs for the 500‐ and 5000‐year return periods of Massawa by considering different GMPEs. The hazard levels estimated by the combined GMPEs tend to lie near the middle of the range defined by the individual GMPEs, and, given the epistemic uncertainty over which model is most applicable, we consider that the combined model is the best choice for the hazard calculations.

We compared our results with previous studies in the region. In general, the studies show similar trends, such as high hazard levels on the plate borders. Obviously, more details can be seen in the present study when compared with Gouin (1976) or Midzi *et al.* (1999). Hagos *et al.* (2006) considered different models; some of which are consistent with our results. However, the PGA values reported in Hagos *et al.* (2006) for the worst‐case scenario are much higher than our predictions. This is mainly due to the choice of the GMPE and the smoothing methods used in addition to the catalog duration.

## Conclusions

We performed a PSHA of Eritrea and the surrounding area using earthquake data compiled from different sources that include the ISC, ESARSWG, and local catalogs. The methods we used here differ from those used in previous comparable studies of the region. The main differences include (1) the data we used cover a longer time span, (2) we used orthogonal regression for the conversion of the magnitudes to a unified scale, which provided a better fit than the standard linear regression, (3) we used adaptive space–time smoothing to construct a seismicity model, (4) we used more recent GMPEs, and a combination of four recent relations with equal weights, (5) we conducted site‐specific hazard assessments for selected important sites within Eritrea, (6) we performed a sensitivity analysis of four input parameters on their influence on hazard levels in Massawa, and (7) we disaggregated the results to determine the magnitude and distance range of seismic sources that dominate the hazard.

The seismic‐hazard maps for return periods of 500, 2500, and 5000 years showed greatest hazard levels around the Red Sea rift, the Gulf of Aden, the Afar depression, and along the western and northern boundaries of the Danakil microplate.

The site‐specific seismic‐hazard analysis of eight selected sites within Eritrea showed a maximum hazard level (at the vibration period of 0.2 s) in Massawa and Bada (0.16*g*) and a minimum hazard level (0.06*g*) in Nakfa at the 500‐year return period. The disaggregation of the hazard for Massawa showed that most of the hazard is contributed by earthquakes of magnitudes 4.5–6 at distances less than 100 km. The sensitivity analysis shows that the most sensitive parameter was the choice of GMPEs and their weight: the Chiou and Young (2008) model yields the highest hazard levels, whereas the Campbell and Bozorgnia (2008) model yields the lowest hazard levels. We recommend using the four relations with equal weights.

The eastern and southeastern parts of Eritrea are seeing large infrastructure development (e.g., in the Massawa and Bada regions) and these areas lie within the high hazard levels in our maps. The seismic risk associated with such development can be significant and these infrastructures need to be developed responsibly and sustainably. The majority of Eritrean buildings, particularly in the countryside, are constructed with low quality and poorly reinforced masonry. This makes the buildings vulnerable even to moderate‐magnitude earthquakes. The damaging 1921 earthquake around Massawa (Gouin, 1979; Ogubazghi and Woldehaimanot, 1998) and the 2011 earthquake of Siroru (Ogubazghi *et al.*, 2014) have already attested to this. As a result, we believe seismic‐hazard assessment of the country is a crucial prerequisite for sustainable economic development and to protect the population at risk.

We recognize the uncertainties in our results, which largely arise from the lack of homogeneous detection levels for small‐magnitude earthquake catalogs, appropriate GMPEs, and fault geometry information. However, we hope our results will stimulate new research in this area and that the hazard maps can be updated when more information is available. The establishment of a seismic network in Eritrea would also allow better characterization of small‐magnitude earthquakes in the region. With this in mind, we hope that our key findings will inform the development of a new building code for Eritrea and contribute to the sustainable economic development of the country.

## Data and Resources

Data sources for our assessment include that of International Seismological Centre (ISC) and Global Centroid Moment Tensor (CMT). The ISC catalog was searched through www.isc.ac.uk/iscbulletin/search/catalogue/ (last accessed December 2014) and the Global CMT database was searched through http://www.globalcmt.org/CMTsearch.html (last accessed January 2016). We also used data obtained from Eastern and Southern African Regional Seismological Working Group (ESARSWG), Eritrea local database, historical data, and events located from our deployment. These data can be obtained through B. G. Information on the potash mining plan around Bada was taken from http://www.danakali.com.au (last accessed May 2016).

## Acknowledgments

B. G. is funded through a Ph.D. scholarship by the University of Bristol and Engineering and Physical Sciences Research Council (EPSRC; Grant Number DTG EP/L504919/1). We would like to thank Eastern and Southern African Regional Seismological Working Group (ESARSWG) for providing us with their catalog that we used in our analysis. The local earthquake catalog during 2011–2012 was generated from data collected using the facilities of Seismic Equipment Infra‐Structure in the UK (SEIS‐UK) supported by the Natural Environment Research Council under Agreement R8/H10/64. Fieldwork support was provided by the Natural Environment Research Council (NERC; Award Number NE/J012297/1 “Mechanisms and implications of the 2011 eruption of Nabro volcano, Eritrea”). C. O. acknowledges support from the NERC Centre for the Observation and Modelling of Earthquakes, Volcanoes and Tectonics. D. K. is supported by NERC Grant Number NE/L013932/1 and Grant Number OSR‐2015‐CRG4‐2643 from King Abdullah University of Science and Technology; and F. I. K. is funded through NERC studentship NE/L002531/1 and a grant to Graduate School of the National Oceanography Centre of Southampton (GSNOCS) from Roy Franklin O.B.E. We thank the two anonymous reviewers and Associate Editor Mark Stirling for their valuable comments on the original article.

- Manuscript received 2 July 2016.