Updated seismic hazard curves, maps, and spectra for the northern Dominican Republic using a probabilistic seismic hazard analysis

This article presents updated seismic hazard curves, spectra, and maps of ground motion intensity measures for the northern region of the Dominican Republic (DR) obtained using a probabilistic seismic hazard analysis (PSHA). The analysis performed uses as input data an earthquake recurrence model based on fault slip rates derived from GPS measurements published in the aftermath of the 2010 Haiti earthquake. The seismicity rate data are used to calibrate a composite characteristic earthquake model, which is combined with a Poisson process to provide a temporal characterization of earthquake occurrence. The seismic hazard curves and maps presented include parameters such as (horizontal) peak ground acceleration and pseudo-spectral response accelerations at 0.2s and 1.0s periods for 5% damping at firm rock sites. The results show that the ground motion parameters with a 2% probability of exceedance (PE) in 50 years determined in this study are up to 46% larger than the corresponding parameters specified in the current DR building code seismic hazard maps for the northern DR. Moreover, the design response spectra for a site in the city of Santiago specified in the code is significantly lower than the 2% PE in 50 years uniform hazard spectra determined in this study for vibration periods smaller than 0.5s, a range that includes the majority of the structures that define the built environment of the DR.


Introduction
The Hispaniola Island is located in a seismically active region bounded by two thrust systems defined by the convergence of the North American and Caribbean plates. Two strike-slip faults, the Septentrional fault and the Enriquillo-Plantain Garden fault, bound the internal Gonâve microplate as shown in Figure 1. The potential of the Hispaniola's seismic faults to cause destructive earthquakes that result in severe damage and losses was demonstrated by the 2010 and 2021 Haiti earthquakes, the former considered one of the most devastating disasters caused by a natural event in recent history (Cavallo et al. 2013). The 2010 Haiti earthquake (January 12, 2010) had a 7.0 M W magnitude, while the 2021 Haiti earthquake (August 14, 2021) had a 7.2 M W magnitude. The 2010 Haiti earthquake caused a coastal uplift with no significant surface rupture along the Enriquillo-Plantain Garden fault (Hayes et al. 2010). The earthquake resulted in major life and economic losses with a death toll estimated to be 300,000 casualties and an economic loss estimated to be nearly US$14 billion, exceeding the gross domestic product of the country DesRoches et al. 2011).
In addition to the seismic activity along the Enriquillo-Plantain Garden fault in the southwestern region of the Hispaniola, the northern region of the island has also long been recognized as a high seismic hazard zone due to the Septentrional fault (Calais et al. 1998;DeMets et al. 2000;Calais et al. 2002;Mann et al. 2002). The Septentrional fault is a left-lateral strike-slip fault that crosses the northern region of the Hispaniola and accommodates part of the motion between the North America and Caribbean plates; a segment of the fault runs close to several densely populated cities of the Dominican Republic (DR). Trenches excavated in the area and radiocarbon analyses have shown that the fault has accumulated at least 5m of deformation and that the most recent severe surface rupture in the central portion of the fault likely occurred before 1260 (Mann et al. 1998). The recurrence interval of the fault is estimated to be between 800 and 1200 years (Prentice et al. 2003), and the slip rate is estimated to be 12 mm/year with a characteristic magnitude or maximum moment magnitude of rupture of 7.8 Frankel et al. 2011). Based on GPS and paleoseismological slip rate estimates, the Septentrional fault may be in the late phase of its rupture cycle .
The Septentrional fault is the main contributor to the seismic hazard and risk in the northern region of the DR due to its proximity to densely populated cities such as Santiago and Puerto Plata (see Figure 2). It is of interest to study and quantify the characteristics of ground motions induced by a potential rupture of this and other nearby faults, and its effects on the built environment (Erazo 2019;Erazo 2020;Rojas-Mercedes et al. 2020;Erazo and Taveras 2021).
The quantitative characterization of earthquake effects at a site or a region is known as a seismic hazard analysis (SHA). The first SHA proposed in the literature was based on a deterministic analysis performed using an estimate of the maximum earthquake that a source can generate (traditionally referred to as the maximum credible earthquake) and the closest distance from a site of interest to the source. Some of the limitations of a deterministic SHA based on a "worst-case scenario" include the time-independent nature of the estimated parameters and the difficulty of rigorously combining the effects from multiple sources considering the time dependency (recurrence) of ground motions generated by each source.
The use of a probabilistic seismic hazard analysis (PSHA) was proposed by Cornell (1968) to address some of the aforementioned limitations of deterministic seismic hazard analyses. A PSHA provides a rigorous and consistent framework to estimate earthquake-induced ground motion parameters at a site, considering the hazard contribution from multiple sources accounting for the likelihood of each individual source to generate earthquakes at different rates. The objective of PSHA is to determine the probability (or rate) of exceeding various ground motion intensity measures, considering the main sources of uncertainty identified by the analyst (McGuire 1995;Rathje and Saygili 2008;Han and Davidson 2012;Rodriguez-Marek et al. 2014;Abrahamson et al. 2019;Gerstenbserger et al. 2020;Baker et al. 2021). PSHA allows the direct treatment of uncertainties such as the likelihood of earthquake magnitudes that can be generated by a source, as well as uncertainties related to rupture length, distance from a site to potential rupture locations, and ground motion prediction equations, among others. Despite some of the challenges involved in defining the inputs to probabilistic seismic hazard analyses (such as limited strong earthquakes data for some regions), the outputs of PSHA have shown to be robust and useful in engineering practice when applied to the design of civil structures, helping to protect public safety against earthquakes. PSHA is the most widely accepted tool to quantify earthquake ground motion intensity measures and its predictions have been tested with observations (Kammerer and Ake 2012).
This article presents the results of a PSHA performed for the northern Dominican Republic, a region considered of high seismic hazard and for which PSHA outputs such as hazard curves and uniform hazard spectra are lacking in the literature. The analysis accounts for the contribution to the seismic hazard of the main sources of known significant seismicity within a 200-km range from the sites of interest. The seismic sources considered in this study are as follows (refer to Figure 2): the Septentrional fault-(3)SFZ, the Northern Hispaniola Subduction Zone-(1)NHFZ, the Matheux-Neiba fault-(10)NFZ, the Enriquillo-Plantain Garden fault-(11)EPGEZ, and the Muertos Trough Subduction Zone-(14)MTFZ. No background seismicity is considered, and although excluding background seismicity underestimates the overall seismic hazard, the impact is typically negligible at the hazard level of a 2% probability of exceedance in 50 years (2475-year return period) which is the focus of this study. To the best knowledge of the author, this is the first peer-reviewed PSHA study in the literature that includes PSHA outputs specifically for densely populated cities in the northern DR such as Santiago and Puerto Plata.
The PSHA presented herein uses as input data fault slip rates derived from GPS measurements for the region published in the aftermath of the 2010 Haiti earthquake (Frankel et al. 2011). The rates are used to calibrate a composite characteristic earthquake model for earthquake recurrence of the type proposed by Youngs and Coppersmith (1985), in which the neighborhood of the characteristic earthquake magnitude is modeled as a uniform random variable and the low-moderate magnitude range as an exponential tail. The model is combined with a Poisson time-independent model for earthquake occurrence. No time correlation or hazard from aftershocks and foreshocks is considered.
The ground motion intensity measures of interest include horizontal peak ground acceleration (PGA), and pseudo-spectral (absolute) response accelerations (PSA) at 0.2s and 1.0s periods for 5% damping at firm rock (class B) sites. The rest of the paper is organized as follows. The next sections discuss the recurrence model and ground motion prediction models employed in this study for spatio-temporal characterization of earthquakes and the resulting ground motions. Then, seismic hazard curves are presented for two sites located in Santiago and Puerto Plata, including the uniform hazard spectra with a 2% and 10% PE in 50 years for both sites. Finally, seismic hazard maps with a 2% PE in 50 years (2475-year return period) for the northern region of the DR are presented for PGA and PSA parameters; the maps are useful to define the design spectra used in engineering practice for the structural design of civil infrastructure. The results are compared to the current design parameters specified in the DR building code. The PSHA results presented can be employed to update outdated DR seismic hazard maps and as a baseline for cross-correlation with future seismic hazard studies for the region, particularly the northern DR.

Magnitude-recurrence model
The recurrence of earthquakes in a region can be quantified using the mean cumulative annual rate of exceedance (also referred to as the cumulative rate of earthquakes or frequency of earthquakes) λ m defined as the ratio of the number of earthquakes of Moment magnitude (M W ) greater than m, divided by the time interval under consideration. The Gutenberg-Richter (GR) law was proposed in 1944 to model the magnitude-recurrence relationship of earthquakes in a region of interest, resulting in an exponential probability distribution for earthquake magnitude.
The GR model is unbounded with no maximum earthquake magnitude that can be generated at a fault. In order to address this drawback, the model can be modified to include lower and upper bounds to limit the earthquake magnitude to an interval, resulting in a truncated exponential probability model for earthquake magnitude. A limitation of the bounded GR model is that it tends to underestimate the observed rate of occurrence of large earthquakes and geological data (Youngs and Coppersmith 1985). In particular, paleoseismicity studies suggest that individual faults periodically generate earthquakes of approximately the same magnitude (characteristic earthquake) and that this magnitude is the maximum that can be generated by the fault; the GR model tends to underestimate the rates near the characteristic magnitude for faults with a characteristic earthquake (Youngs and Coppersmith 1985). The composite characteristic earthquake model proposed by Youngs and Coppersmith attempts to address this limitation by modeling the low-moderate magnitude range with an exponential tail and the magnitude range close to the characteristic earthquake as a uniform (constant) distribution. The transition between the exponential and uniform regions takes place at a magnitude m c = m char − Δm c with the exceedance rate of the uniform region defined as the rate at m ′ = m c − Δm ′ where m char is the fault characteristic magnitude, and Δm c and Δm ′ are model parameters (Youngs and Coppersmith 1985;Gupta 2007).
Given that earthquake recurrence rates vary among seismic sources a recurrence model needs to be calibrated for each source. This study considers the following seismic sources for the PSHA: the Septentrional fault, the Northern Hispaniola Subduction Zone, the Matheux-Neiba fault, the Enriquillo-Plantain Garden fault, and the Muertos Trough Subduction Zone. The characteristic earthquake magnitude, slip rates, and predicted earthquake annual exceedance rates for the sources obtained by Frankel et al. (2011) based on GPS data constrained by historic large earthquakes are shown in Table 1; this data is based on studies performed in the aftermath of the devastating 2010 Haiti earthquake, and it is used herein as input to calibrate a recurrence model for the sources. It is worth pointing out that the resulting recurrence models are based on estimated slip rates and seismic  (Frankel et al. 2011). Moreover, no earthquake of magnitude greater than M w 6 has occurred in the DR since the 2003 Puerto Plata earthquake, which implies that the slip rates used are the most updated ones available in the literature. The earthquake magnitude recurrence rates that result from using the (incomplete) instrumental seismicity data available for the DR are significantly higher than those presented in Table 1 (see Appendix) and are inconsistent with the rates reported in the literature from paleoseismicity studies (Frankel et al. 2011). It will be shown in the following section that the seismic hazard in the northern region of the DR is dominated by the Septentrional fault due to a combination of high large magnitude recurrence rates and its close distance to sites in the region; for this reason, the recurrence model of the Septentrional fault will be discussed with more detail. The exceedance annual rates for the Septentrional fault for magnitudes of 6.5 and 7 are, respectively, λ 6.5 = 0.021 and λ 7.0 = 0.0078 (Frankel et al., 2011). The rates are based on a model where 67% of the seismic moment is released by characteristic earthquakes, and the remaining 33% of the seismic moment follows an exponential tail. The aforementioned exceedance rates and the fault slip rate (12mm/year) are used as constraints to calibrate a characteristic earthquake model of the composite type proposed by Youngs and Coppersmith (1985) for the fault. The model calibration was performed using an optimization procedure constrained by the recurrence rates and slip rates, with tuning parameters Δm c and Δm ′ ; the calibration best-fit parameters were Δm c = 0.35 and Δm ′ = 0.35. The resulting model is shown in Figure 3, which shows the model recurrence rates and the corresponding return period. The minimum rate considered in the model is λ min = λ 6.5 and the maximum rate is the characteristic earthquake rate. The adopted characteristic earthquake model is compared to the Gutenberg-Richter model employed in the current DR building code in the Appendix.
Using the same procedure described before a composite characteristic model was calibrated for the Matheux-Neiba fault and the Enriquillo-Plantain Garden fault. Given the lack of additional information, the model adopted for the Northern Hispaniola and Muertos Trough Subduction Zones consisted of a characteristic-only model with the rates defined in

Seismic hazard characterization in a time interval
The magnitude-recurrence model discussed in the previous section allows us to estimate the recurrence of earthquakes and the earthquake magnitude probability distribution at a site without considering a time interval or exposure time in the analysis. For engineering applications, it is also of interest to characterize the potential occurrence of earthquakes in a given time interval that involves the lifetime of a civil structure or an exposure time (e.g., 50 years, 100 years, etc.). After the exposure time is defined, the objective is to determine the probability of exceeding an earthquake magnitude or ground motion intensity measure such as peak ground or response spectral accelerations in order to design civil structures for a parameter value with a target "small" probability of being exceeded in the defined time interval, or equivalently, with a target return period. Defining an exposure time implies that a model for earthquake occurrence in time is needed. The Poisson process is typically adopted to model events related to earthquake occurrence in time. In this model, the number of occurrences X T of an event in the time interval [0,T] is modeled as a Poisson distribution. The variable X T can be, for example, the number of earthquakes that exceed a given magnitude, or the number of times that a ground motion intensity measure exceeds a given value. The probability of exceedance (PE) is defined as the probability that X T occurs at least once in the time interval considered, and it is given by PE ≡ p(X T ≥ 1) = 1 − e −λT where λ is the mean annual rate of occurrence of X T .
It is worth pointing out that modeling earthquake occurrence as a Poisson process is not in agreement with some of the physical processes that generate earthquakes and the elastic rebound theory, since it assumes that earthquakes generated by a source are stochastically independent. Given the strain energy accumulation and release mechanisms that take place at a fault, there is a time correlation (clustering) among earthquakes, which implies that the mean occurrence rate is time-dependent (Ogata 1999; Cramer et al. 2000;Sykes and Menke 2006;HERP 2010). Despite this limitation of the Poisson process model, the results of PSHA have shown to be robust to the recurrence model (including timedependent models) and provide a systematic framework to define engineering design parameters, which justifies its use in earthquake engineering applications (Chang et al. 2017).

Ground motion intensity measures and models
A fundamental objective of PSHA is the estimation of ground motion intensity measures at a site and their associated recurrence. The ground motion characteristics that have the most significant impact on the response and behavior of civil structures are duration, amplitude (intensity), and frequency content. Some of the factors affecting these parameters are as follows: Type of seismic fault and fault orientation, rupture depth, and rupture direction with respect to the site (directivity).
Earthquake magnitude, which in turn depends on the area of the ruptured zone, fault displacement, and shear modulus of the rock adjacent to the fault, among other parameters.
Characteristics of the propagation medium include local soil properties (such as wave velocity propagation) and topography.
Given the difficulty of analytically modeling the ground motions induced by earthquakes using mechanics-based models, empirical models known as attenuation models or ground motion prediction equations (GMPE) are typically used to perform this task. GMPE are models developed using regression analyses and past seismicity data to relate predictor variables, such as earthquake magnitude and distance, to ground motion intensity measures . When sufficient seismicity data is available GMPEs can be calibrated for a specific site of interest. In some instances, however, this is not the case particularly in regions where instrumental seismicity data is scarce given the large recurrence rate of large earthquakes, and thus data from different regions with similar fault mechanisms and local characteristics need to be combined.
In this study, two sets of GMPE are used for crustal faults and subduction zones. For crustal earthquakes the NGA-West2 GMPE is adopted ). The models were developed using measured data from an expansion of the PEER (Pacific Earthquake Engineering Research Center) strong motion database recorded from shallow crustal earthquakes in active tectonic regions. The predictor variables include magnitude scaling, geometric attenuation, type of fault, hanging-wall effects, shallow linear/nonlinear site response, basin response, hypocentral depth, and rupture dip. The models used (with equal weights) are the models by Abrahamson et al. (2014)  For earthquakes on the interface of subduction zones, the GMPE employed are the models by Zhao et al. (2006), Atkinson and Boore (2003), and Youngs et al. (1997) weighted by factors of 0.5, 0.25, and 0.25, respectively; a higher weight is typically assigned to the Zhao et al. model given the larger database used in its calibration, particularly for distances smaller than 100km. The weights are based on the 2008 USGS weights used for the U.S. national seismic hazard maps (Petersen et al. 2008). It is worth stressing that the DR does not have an updated database of measured strong ground motions in order to validate or calibrate GMPEs specifically for the northern DR. Despite this drawback, the results presented in the following sections showed to be robust to variations in the weight distribution assigned to the adopted GMPEs.
GMPEs are parameterized by a site-to-rupture zone distance that accounts for radiation damping or seismic wave attenuation. For nearly vertical faults (such as the Septentrional fault), the Joyner-Boore distance (closest distance from the surface projection of the rupture zone to the site) is generally adopted, and defined as r JB = √ r 2 e + r 2 min where r min is the smallest distance from the site to the surface projection of the fault, and r e is the epicenter location measured with respect to the point in the surface projection of the fault closest to the site. The epicenter location r e is modeled as a uniform random variable along the length of the fault, implying that the likelihood of a rupture initiation at different locations along the fault is constant. This is equivalent to floating or moving the rupture zone assuming that the likelihood of the epicenter location is constant along the longitudinal direction of the fault. For GMPE that uses the closest distance to the rupture plane, r rup , the distance metric is given by (for nearly vertical faults) r rup = √ r 2 jb + Z 2 TOR where Z TOR is the depth to the top of rupture; a probability distribution for Z TOR is defined as a function of magnitude using the results by Wells and Coppersmith (1994) and Kaklamanos et al. (2011). For dipping fault planes in subduction zones, the closest distance of the rupture to the site is employed.
The GMPEs are used for global regions and uniform site conditions with an average shear-wave velocity v s30 of 760 m/s, which corresponds to Class B firm rock sites; the rest of the GMPE parameters are selected as the recommended values defined for each model (Gregor et al. 2014). The analysis requires the characteristics of the fault rupture, such as depth to the top of the rupture area, hypocenter depth, and rupture length; parameters related to fault rupture are estimated using the results by Wells and Coppersmith (1994) and Kaklamanos et al. (2011).

Seismic hazard curves for ground motion intensity measures
The main output of a PSHA is the rate (or probability) of exceeding different ground motion intensity measures at a site or region. The parameters of interest in engineering applications typically include (horizontal) peak ground acceleration and pseudo-spectral response accelerations for different vibration periods and 5% damping. In particular, the spectral accelerations at 0.2s and 1s have traditionally been used in building codes to define the design response spectra and seismic loads for the design of civil structures (Villaverde 2009). For faults modeled as a line source, the annual rate of exceedance of the event Y > y (related to exceeding various ground motion intensity measures) is given by where P(Y > y | m, r) is the conditional probability of exceeding the ground motion level y given the occurrence of an earthquake of magnitude m at a distance r and it is defined by the GMPE, P(m) is the probability density function (PDF) of earthquake magnitude obtained from the recurrence model (see Appendix), and P(r | m) is the PDF of the seismogenic site-to-rupture location distance conditional on earthquake magnitude; as mentioned before, the seismogenic distance r and the earthquake magnitude m (1) are related and thus stochastically correlated random variables.
The rate λ min is the rate of earthquake magnitudes greater than the minimum magnitude considered in the analysis. When multiple sources are within the range of a site of interest, the total rate is computed by adding the rates estimated for the individual sources. This procedure allows to rigorously combine the parameters/events rates of different sources. The resulting annual rate of exceedance λ Y > y is used to determine the probability of occurrence of the event Y > y in a time interval or exposure time.

Results for Santiago and Puerto Plata
The seismic hazard curves for different ground motion intensity measures are presented next for sites near the center of the two main cities in the northern DR, Santiago and Puerto Plata. Based on their population and built environment, these are the cities with the highest seismic risk in the country. Figure 4 shows the seismic hazard curves for (horizontal) peak ground acceleration (PGA) for a site close to the center of Santiago (Figure 4); the site is located at a distance r min = 8 km from the Septentrional fault, measured to the surface projection of the fault. The hazard curves show the contribution of all the sources to the total hazard; as can be seen, the total hazard curve is dominated by the Septentrional fault. This is attributed to the close distance between the Septentrional fault and Santiago, with the second closest source being the Northern Hispaniola Subduction Zone located at a distance r min = 62 km; moreover, the latter source has a significantly lower recurrence rate (see Table 1). Figure 4 also shows the sensitivity of the hazard curve based on the Septentrional fault to the four GMPE employed in the PSHA (with equal weights) for crustal faults.

Results for a site in Santiago
The PGA values with a 2% and 10% PE in 50 years are, respectively, 0.83g and 0.50g. According to the DR building code, the corresponding PGA values with 2% and 10% exceedance probabilities for this site are 0.60g and 0.43g (MOPC 2011). The PGA corresponding to a 2% PE in 50 years obtained in this study is 38% higher than the value specified in the DR building code. The 2% PE in 50 years PGA value estimated in this study is in agreement with the results in Frankel et al. (2011), which estimate the corresponding PGA for the site to be in the range of 0.80-1.00g.
As mentioned before, the pseudo-spectral (absolute) acceleration (PSA) for different vibration periods and 5% damping are used to define the design spectra in building codes. The design spectrum is traditionally defined using the PSA values at periods of 0.2s (PSA 0.2 ) and 1s (PSA 1 ) for class B firm rock sites, and then it is further modified by a set of factors   to account for local site effects such as soil amplification. In some building codes (such as ASCE 07), the spectrum with a return period of 2475 years is further reduced by a 2/3 factor to account for overstrength and a life safety performance limit state. For consistency in the results, the PSHA spectral ordinates determined in this study will not be reduced by the 2/3 factor throughout the manuscript, and thus, the 2475-year PSHA spectra from this study can be directly compared to the 2475-year spectra specified in the DR building code. Figure 5 shows the seismic hazard curves for PSA 0.2 for the same site in the center of Santiago. The PSA 0.2 value with a 2% PE in 50 years is 2.05g, while the PSA 0.2 value with a 10% PE in 50 years is 1.15g. According to the DR building code, the PSA 0.2 values with 2% and 10% probability of exceedance for this site are, respectively, 1.40g and 1.00g (MOPC 2011). The PSA 0.2 corresponding to a 2% PE in 50 years obtained in this study is 46% higher than the value specified in the DR building code. The 2% PE in 50 years PSA 0.2 value estimated in this study is in agreement with the results in Frankel et al. (2011), which estimate the corresponding PSA 0.2 for this site to be in the range of 1.50-2.40g. Figure 6 shows the seismic hazard curves for PSA 1 . Based on the analysis, the PSA 1 value with a 2% PE in 50 years is 0.69 g; similarly, the PSA 1 value with a 10% PE in 50 years is 0.36g. According to the DR building code, the PSA 1 values with 2% and 10% exceedance probabilities for the site are, respectively, 0.72g and 0.48g (MOPC 2011). The results are summarized in Table 2.
Using the procedure described before, the seismic hazard curves for PSA at different vibration periods can be readily obtained. A uniform hazard spectrum (UHS) is formulated by combining the PSA values at different vibration periods with the same return period or probability of exceedance in an interval (Loh et al., 1994;Wen and Wu, 2001;Gavin and Dickinson, 2011;Lee et al., 2013) The UHS with a 2% PE in 50 years is used to define the seismic loads employed for the seismic design of most conventional civil structures. This design spectrum is referred to as the maximum considered earthquake (MCE) in building codes, although the terminology is misleading given that the output of the PSHA is an envelope or convolution that considers the likelihood of different earthquake magnitudes and the potential ground motion characteristics induced by these earthquakes as estimated by the GMPE. Figure 7 shows the uniform hazard spectrum for a 2% PE in 50 years determined for the site under consideration in the city of Santiago. The figure also shows the 2% PE in 50-year spectra for the site based on the specifications in the current seismic building code of the DR for rock sites; both spectra correspond to class B firm rock sites with no  additional factors applied in order to compare the spectra with a 2475-year return period directly. As can be seen, the UHS obtained in this study significantly exceeds considerably the spectra specified in the DR building code for vibration periods smaller than 0.5s. It is worth noting that this is the range of the vibration periods of most low-rise structures and buildings under five stories that constitute the majority of the built environment of the city of Santiago and the rest of the country. The results imply that these structures are currently being designed for significantly smaller forces than the design forces resulting from a response spectrum with a 2% PE in 50 years. As can be seen, the low-frequency (high period) range of the spectra is consistent. This is attributed to the fact that the Campbell 1997 model employed in the current DR code predicts significantly larger spectral accelerations than the NGA West2 models for vibration periods larger than 0.5s, compensating for the overall reduction in seismic hazard that results from the recurrence model as further discussed in the Appendix. For vibration periods smaller than 0.5s, the median parameters estimated by the Campbell 1997 (C97) and the NGA West2 models are consistent, but the larger standard deviation associated with the NGA West2 models, combined with the adopted characteristic earthquake recurrence model, results in an increase in the seismic hazard in this region of the spectra.    The difference in the results presented herein and the parameters specified in the DR building code is attributed to two main factors: The spectral parameters in the DR building code are based on an analysis that employs an unbounded GR model for earthquake recurrence, which significantly underestimates the probability of occurrence of large magnitude earthquakes associated with large PSA values (see Appendix). Moreover, the recurrence model was calibrated using only past seismicity data which in addition to being scarce, included mostly data of earthquake magnitudes in the 3-5 M W range (SODOSISMICA, 2009). The resulting probability density function shows that earthquakes of magnitude M W > 6 (and the associated ground motion parameters induced) receive a small weight when integrated with Equation 1 since they have a small probability associated with them. The characteristic earthquake model employed in this manuscript is based on earthquake rates that include information related to the fault slip rate, as well as historic records used to constrain the occurrence rates of strong earthquakes which results in a significantly higher probability associated with large earthquakes. We believe that our model captures the seismicity of the fault with more accuracy than the unbounded GR model based on small amplitude waveforms.
The parameters in the DR building code are based on the ground motion prediction equation proposed by Campbell in 1997(Campbell, 1997, which is outdated and has been significantly refined over the last 20 years using strong motion data from a significantly larger earthquakes database . The standard deviation of the ground motion intensity measures estimated by recent models is considerably larger than the estimates provided by the Campbell 1997 model (see Appendix). The smaller uncertainty associated with parameters predicted by the Campbell 1997 model implies that ground motion parameter values larger than the median (which are associated with strong shaking) receive a small weight when integrated into Equation 1, resulting in a small probability and thus a small rate for such parameter values.
These factors contribute to a decrease in the estimates of ground motion parameters associated with large earthquakes. The UHS in Figure 7 stresses the importance to update the current seismic hazard specifications in the DR building code, especially for the northern region of the country. It is worth pointing out that an analysis with uniform hazard as a target does not result in a uniform probability of collapse (uniform risk) at different locations as shown in Li et al. (2010) for light-frame wood residential structures due to uncertainty in the capacity (strength) of the structure and variability in the shape of the hazard curves. To achieve a uniform risk-based design, the parameters are typically adjusted by a factor between 0.85 and 1.15 (Luco et al., 2007). Figure 8 shows the seismic hazard curves for (horizontal) peak ground acceleration (PGA) for a site close to the center of Puerto Plata. The site is located at a distance r min = 26 km to the Septentrional fault, and a distance r min = 26 km to the Northern Hispaniola Subduction Zone. The PGA value with a 2% PE in 50 years for the site is 0.50g. Similarly, the PGA with a 10% PE in 50 years is 0.30g. The curves show the contribution of the various sources to the total hazard, and as can be seen, the Septentrional fault and the Northern Hispaniola Subduction Zone have a comparable contribution to the total hazard in Puerto Plata. The PGA for various hazard levels in Puerto Plata is not specified in the DR building code.

Results for a site in Puerto Plata
The seismic hazard curves for PSA 0.2 are shown in Figure 9. The PSA 0.2 value with a 2% PE in 50 years is 1.20 g. Similarly, the PSA 0.2 with a 10% PE in 50 years is 0.70g. The PSA hazard values are also not Instead, the code suggests using the values for Santiago, assuming that these parameters are higher than the corresponding Puerto Plata parameters and thus provide an upper bound. Figure 10 shows the seismic hazard curves for PSA 1 . Based on the analysis results, the PSA 1 value with a 2% PE in 50 years is 0.43g, while the PSA 1 value with a 10% PE in 50 years is 0.22g. The PSA 1 hazard values are also not specified in the code. The results are summarized in Table 2. Figure 11 shows the uniform hazard spectra with a 2% PE in 50 years determined in this study for a site in Puerto Plata. As mentioned before, a 2% PE in 50 years design spectrum is not defined in the DR building code for this site, and the code recommends using the design spectrum for Santiago. It is worth pointing out that in 2003 Puerto Plata was struck by a 6.4 M W earthquake that caused significant structural damage and the collapse of multiple houses and lowrise buildings (Dolan and Bowman, 2004). In particular, several public schools showed partial or total collapse; fortunately, the earthquake occurred afterhours and the school buildings were empty.

Seismic hazard maps
The uniform hazard spectra for a site in the cities of Santiago and Puerto Plata were presented in the previous section. Given that a site of interest can be in any spatial location within a region, hazard maps that provide ground motion measures with a constant return period (or probability of exceedance in a period) provide a useful tool for implementation in building codes. The maps are defined based  on contours or lines of constant parameter values, which makes them suitable for engineering practice. Figure 12a shows the hazard map for PGA with a 2% PE in 50 years and firm rock sites based on the PSHA performed in this study, while Fig. 12b shows the corresponding PGA hazard map specified in the current DR building code. As can be seen, the hazard map in the building code shows significantly lower Fig. 12 Seismic hazard map for PGA with a 2% probability of exceedance in 50 years; a this study, and b DR building code (MOPC 2011) hazard levels than this study, especially in the nearfault regions. It is worth pointing out that based on the results of Frankel et al. (2011) the PGA value with a 2% PE in 50 years along the Septentrional fault is in the 1.0-1.8-g range; the PGA value of 1.15g determined in this study is consistent with this range, Fig. 13 Seismic hazard map for PSA 0.20 with a 2% probability of exceedance in 50 years; a this study, and b DR building code (MOPC 2011) while the 0.64g value specified in the code is significantly smaller than range lower bound. Figure 13 and Figure 14 show, respectively, the hazard maps for PSA 0.2 and PSA 1 with a 2% PE in models used to develop these maps. According to the results of Frankel et al. (2011), the value of the PSA 0.2 parameter along the Septentrional fault is in the 2.4-3.2-g range, which is in agreement with the 2.85-g value determined in this study. The corresponding value specified in the DR building code is 1.51g, which again is significantly smaller than the range lower bound. Note that the seismic hazard maps presented do not have a deterministic upper limit near the faults, since this is typically defined when adopted in building codes. For example, the ASCE 7-22 maps use a factor of 1.8 times the estimated median response to the characteristic earthquake of the faults to define the upper limit (ASCE, 2022).
The results of this study could potentially be used by regulating authorities of the Dominican Republic as a baseline to update the current seismic hazard maps specified in the DR building code, Moreover, the results provide a rigorous basis for comparison of future seismic hazard studies in the Dominican Republic.

Conclusions
This paper presented the results of a probabilistic seismic hazard analysis (PSHA) performed for the northern region of the Dominican Republic (DR). The PSHA uses as input recurrence models calibrated using fault slip rates derived from GPS measurements published in the aftermath of the 2010 Haiti earthquake. The recurrence models were combined with a Poisson process model and a set of recently proposed ground motion prediction equations (GMPE) for spatio-temporal characterization of ground motion intensity measures at a site. The PSHA performed accounted for uncertainties related to earthquake magnitude, distance from a site to the potential rupture locations along the fault, and the aleatory uncertainty associated with the regression analyses performed to obtain the GMPE.
The estimated ground motion parameters included horizontal peak ground acceleration (PGA) and pseudo-spectral (absolute) response accelerations at 0.2-s (PSA 0.20 ) and 1.0-s (PSA 1 ) periods for 5% damping in firm rock (class B) sites. A set of seismic hazard curves were developed for the PGA, PSA 0.20 , and PSA 1 parameters for two sites in the center of the cities of Santiago and Puerto Plata. The results included the hazard spectra with a 2% probability of exceedance (PE) in 50 years for both sites. The results showed that the ground motion parameters such as PGA and PSA for vibration periods smaller than 0.5s defined in the current DR building code are inconsistent with a 2% PE in 50 years or 2475-year return period hazard level. This implies that low-rise civil structures with vibrations periods smaller than 0.5s are currently being designed for seismic forces significantly smaller than the target forces expected from a design response spectrum with a 2475-year return period. Structures in this class include one to five stories houses and buildings, which represent the majority of the built environment of the DR.
In addition to hazard curves and spectra, hazard maps with a 2475-year return period were developed for the PGA, PSA 0.20 , and PSA 1 parameters for the northern region of the Dominican Republic. The maps were compared to the current hazard maps specified in the DR building code. It was shown that the response parameters specified in the current code are underestimated and inconsistent with the 2475-year hazard level from this study throughout the northern region of the DR. The results presented herein showed to be consistent with the results by Frankel et al. (2011) developed for Haiti, but included new PSHA outputs such as uniform hazard spectra and hazard curves for cities in the Dominican Republic. The results highlight the need to update the current design seismic hazard maps of this country given the high strain accumulation levels in the Septentrional fault reported in the literature.
Future work is aimed at including spatially smoothed (background) seismicity to account for earthquakes not attributed to crustal faults and subduction zones, and extending the analysis to the rest of the country. motion prediction equations (GMPE) employed for spatio-temporal characterization of ground motion parameters.
To illustrate this, Figure 15 shows the earthquake recurrence models used in this study (a composite characteristic earthquake CE model) and the current DR building code (an unbounded Gutenberg-Richter GR model with parameters a = 1.98 and b = 0.49) for the Septentrional fault (SODOSISMICA, 2009); the figure also depicts the resulting probability density function (PDF) for earthquake magnitude at the fault. Fig. 15a shows the observed catalog seismicity of the DR, according to the non-profit organization SODOSISMICA (SODOSISMICA, 2009); the data in this catalog is based on earthquakes that occurred mostly within the last 30 years, which is considered insufficient for model calibration purposes as evidenced by the significantly small 0.49 b-value. The GR model (2009) adopted in the DR building code and shown in the figure is based on this seismicity data.
As can be seen, the earthquake rates estimated by the unbounded GR model are larger than those estimated by the characteristic model adopted in this study, particularly in the range close to the characteristic magnitude M W 7.8. The unbounded GR model used in the DR building code was calibrated using a catalog data that included 166 earthquakes in the range 3 ≤ M w ≤ 5 and 8 earthquakes with M w > 6 (SODOSISMICA, 2009); the catalog observation data used in SODOSISMICA (2009) is shown in Fig. 15. Note that although the recurrence rates of the GR model are higher than those estimated by the CE model, the probability of earthquakes of magnitude M w > 6 (and the ground motion parameters associated with them) receive a small weight due to the PDF P(m) shown in Fig. 15b when integrated into the hazard equation (Equation 1); indeed, the PDF of earthquake magnitude based on the GR model is almost two orders of magnitude smaller than the PDF based on the composite characteristic model for the range 6.5 ≤ M w ≤ 7.8 due to the wider domain of the former.
In essence, according to the model used in the current DR seismic code, the majority of the seismic moment is released by earthquakes with magnitude in the range 3 < M w < 6, which in turn are associated with small ground motion parameter values resulting in a lower seismic hazard than the model employed in this study. Similarly, the rate of large ground motion parameters is small since the term P(Y > y) is penalized due to the small weight provided by P(m) for the earthquake magnitude more likely to generate large ground motion response parameters. This chiefly results in lower estimates of the recurrence rate of large ground motion intensity measure values and partially (b) (a)  Fig. 15 Comparison of the recurrence models adopted in the current DR building code (unbounded Gutenberg-Richter GR model) and this study (composite characteristic earthquake CE model); a cumulative earthquake rates, and b probability density functions (PDF) justifies the increase in the seismic hazard of this study with respect to the existing building code. The second major PSHA input that considerably differs from that adopted in the current DR building code is the set of GMPE employed. Figure 16 depicts the GMPE median and standard deviation for the models used in this study for the Septentrional fault, as well as the model adopted in the current DR code (Campbell, 1997).
As can be seen, although the median ground motion intensity measures are within the same range, the standard deviation estimated by the NGA-West2 models is considerably larger than the corresponding standard deviation estimated by the Campbell 1997 (C97) model for small distances that are associated with strong shaking. The smaller uncertainty associated with the C97 model implies that ground motion parameter values larger than the median (which are associated with strong shaking) receive a small weight when integrated by P(Y > y | m, r) in Equation 1, resulting in a small probability and rate for such parameter values.
The Campbell 1997 (C97) model is outdated and more recent models (such as the NGA-West2 models) have been developed over the last 20 years using strong motion data from a significantly expanded database. Moreover, only the Campbell 1997 model was used to develop the current DR hazard maps, including regions with subduction zones for which the model does not apply (SODOSISMICA, 2009).

Fig. 16
Comparison of ground motion prediction equations (GMPE) adopted in the current DR building code  and this study for crustal faults (NGA-West2 models); a median PGA, and b PGA standard deviation (logarithmic units)