Risk-informed representative earthquake scenarios for Valparaíso and Viña del Mar, Chile

Different risk management activities, such as land-use planning, preparedness, and emergency response, utilize scenarios of earthquake events. A systematic selection of such scenarios should aim at finding those that are representative of a certain severity, which can be measured by consequences to the exposed assets. For this reason, defining a representative scenario as the most likely one leading to a loss with a specific return period, e.g., the 100-year loss, has been proposed.

We adopt this definition and develop enhanced algorithms for determining such scenarios for multiple return periods. With this approach, we identify representative earthquake scenarios for the return periods of 50, 100, 500, and 1000 years in the Chilean communes of Valparaíso and Viña del Mar, based on a synthetic earthquake catalog of 20 000 scenarios on the subduction zone with a magnitude of M_w≥5.0. We separately consider the residential-building stock and the electrical-power network and identify and compare earthquake scenarios that are representative of these systems. Because the representative earthquake scenarios are defined in terms of the annual loss exceedance rates, they vary in function of the exposed system. The identified representative scenarios for the building stock have epicenters located not further than 30 km from the two communes, with magnitudes ranging between 6.0 and 7.0. The epicenter locations of the earthquake scenarios representative of the electrical-power network are more spread out but not further than 100 km away from the two communes, with magnitudes ranging between 7.0 and 9.0. For risk management activities, we recommend considering the identified scenarios together with historical events.

1 Introduction

Due to the complexity of earthquake events and the response of infrastructure and society to these events, risk managers analyze potential impacts of strong seismic events and test risk management capacities through representative earthquake scenarios (e.g.,  Salgado-Gálvez et al., 2018; Aguirre et al., 2018). Scenario-based analysis enables the modeling and simulation of the complex processes and interactions during and after earthquake events, with a level of detailing that is not possible through a complete probabilistic hazard and risk analysis. As such, the earthquake scenarios are the starting point for such a more detailed risk assessment and for recommendations for improving risk management (e.g.,  Chatelain et al., 1995; Feliciano et al., 2023).

Representative scenarios are commonly selected based on expert knowledge (e.g.,  Aguirre et al., 2018) and past events (e.g.,  Indirli et al., 2011). Synthetic seismic catalogs have also been used for the selection of representative scenarios (McGuire, 1995; Jayaram and Baker, 2009b; Miller and Baker, 2015). A particular approach for scenario selection is based on hazard disaggregation (Bazzurro and Cornell, 1999), which utilizes the conditional probability of different hazard scenarios given an intensity measure (e.g., peak ground acceleration, PGA) at a specific site of interest either equaling or exceeding a threshold (Fox et al., 2016; Fox, 2023). As its name suggests, classic hazard disaggregation does not explicitly consider the losses of the affected engineering systems, which are often a function of the intensity measures at multiple locations and which are subject to uncertainty.

The above concepts were extended to loss disaggregation to find earthquake scenarios in terms of magnitude and hypocentral distance that exceed a loss threshold for residential-building stocks (Goda and Hong, 2009) or infrastructure (Jayaram and Baker, 2009b). Because the spatially accumulated loss can be defined for any portfolio of buildings and infrastructure, loss disaggregation implicitly considers spatially distributed intensity measures. Rosero-Velásquez and Straub (2022) proposed a definition of a representative hazard scenario associated with a loss of return period t, e.g., the 100-year loss, which in general does not correspond to the magnitude or intensity measure of the same return period. It is defined as the most likely scenario that leads to the loss value (i.e., its occurrence) associated with this return period t. They also presented a numerical procedure for selecting the representative hazard scenario in a continuous space of source parameters with a surrogate model and active learning, thus considering the uncertainty in the conditional losses given a hazard scenario. In the context of seismic risk analysis, earthquake scenarios that are representative of t-year loss can be different for different engineering systems, even if they are located in the same area. The definition of Rosero-Velásquez and Straub (2022) differs from the loss disaggregation presented by Goda and Hong (2009) and Jayaram and Baker (2009b) because the latter define the representative scenario as the most likely one to exceed the t-year loss. In this contribution, we compare the two definitions and argue that a definition in terms of the occurrence of the t-year loss is more appropriate for most applications, in line with the findings of Fox et al. (2016) for hazard disaggregation.

Considerable work has been devoted to the study of the seismic hazard, vulnerability, and risk in the Valparaíso coastal area of Chile due to its high population density and economic importance in combination with strong seismic activity. Recent earthquakes that led to significant damage occurred in 1971 with M_w=7.8, in 1985 with M_w=8.0 (Indirli et al., 2011), and in 2010 with M_w=8.8 (de la Llera et al., 2017). Recent studies on the Valparaíso area deal with seismic characterization (e.g.,  Carvajal et al., 2017; Candia et al., 2020), source models (e.g.,  Poulos et al., 2019; Pagani et al., 2021), ground motion models (Montalva et al., 2017), building exposure models (e.g.,  Yepes-Estrada et al., 2017; Jiménez et al., 2018; Gómez-Zapata et al., 2022b, b), damage analysis on individual buildings (e.g.,  Indirli et al., 2011; Jünemann et al., 2015), socio-economic impact (Jiménez Martínez et al., 2020), and the seismic risk analysis of the electrical-power network (Ferrario et al., 2022) and road network (Allen et al., 2022). Additionally, Indirli et al. (2011) identified representative earthquake scenarios using historical events and expert knowledge for generating representative ground motion time series but solely from the hazard point of view and disregarding the risk component.

This paper determines representative earthquake scenarios for different return periods for the residential-building stock and the power supply network in the communes of Valparaíso and Viña del Mar. We adapt and extend the methodology described by Rosero-Velásquez and Straub (2022) for identifying scenarios associated with different return periods from a synthetic earthquake catalog. The representative scenario is found directly by solving a stochastic optimization problem, namely the identification of the mode of the conditional distribution of the source parameters given the occurrence (or exceedance) of the t-year loss among the scenarios in the catalog. The stochastic optimization problem is solved with an active-learning strategy, whereby the uncertainty in the objective function is estimated by bootstrapping.

We introduce the definition of a representative earthquake scenario more formally in Sect. 2. Then we present the methodology for computing the scenarios on a seismic catalog in Sect. 3 and illustrate it with idealized examples in Sect. 4. The description of the study area and the utilized hazard and system models are presented in Sect. 5. The results are given in Sect. 6 and discussed in Sect. 7. In the paper, we employ the notation and acronyms summarized in Appendix A.

2 Definition of a representative earthquake scenario

An earthquake scenario can be described by a vector θ of source parameters, including the magnitude; hypocentral distance; and source longitude, latitude, and depth. In a stochastic model, the scenario is a single realization of a random vector Θ, with joint probability density function (PDF) f_Θ(θ). The PDF of Θ is obtained from one or more seismic source models (e.g., Poulos et al., 2019) and is conditioned on the occurrence of a seismic event, whose frequency (occurrence rate) is λ_H.

An earthquake catalog is a set of n earthquake scenarios of ${\mathit{\theta }}^{\left(\mathrm{1}\right)},\mathrm{\dots },{\mathit{\theta }}^{\left(n\right)}$, which are realizations of Θ. The catalog can be a set of synthetic earthquake scenarios, obtained by random sampling from f_Θ(θ). Alternatively, the catalog can be obtained from past events (e.g.,  Poulos et al., 2019).

Figure 1 Illustration of the representative scenario in the source parameter space, modified after Rosero-Velásquez and Straub (2022), in terms of loss occurrence (a) and in terms of loss exceedance (b).

Synthetic earthquake catalogs have been used in event-based probabilistic seismic hazard analysis (PSHA) and earthquake risk assessment (e.g., Salgado-Gálvez et al., 2018; Ferrario et al., 2022; Allen et al., 2022). The aim of PSHA is to obtain the occurrence rate and distribution of ground motions, taking into account all possible earthquake scenarios (Cornell, 1968; Esteva, 1970), and event-based PSHA utilizes Monte Carlo simulation for sampling earthquake scenarios. Similarly, event-based earthquake risk assessment on spatially distributed systems utilizes synthetic earthquake scenarios for computing the losses, considering the spatial correlation in the ground motion and the vulnerability of the exposed assets (Baker et al., 2021).

For a given earthquake scenario, ground motion models (GMMs) result in spatially distributed intensity measures, e.g., PGA and spectral accelerations, which are input to assess the losses associated with the exposed systems. In the general case, these predictions are stochastic. Thereafter, the model of the engineering system considers the physical and functional vulnerability and results in a loss value L.

Because of the randomness and uncertainty in the earthquake scenario, GMM, vulnerabilities, and exposure, L is a random variable whose cumulative distribution function (CDF) F_L(l) can be obtained by performing an event-based earthquake risk assessment for spatially distributed systems with the synthetic earthquake catalog. By combining this CDF with the earthquake occurrence rate λ_H, one obtains the loss exceedance function λ_L(l):

$$\begin{array}{}\text{(1)}& {\mathit{\lambda }}_{\mathrm{L}}\left(l\right)=(\mathrm{1}-F_{\mathrm{L}}(l\left)\right){\mathit{\lambda }}_{\mathrm{H}}.\end{array}$$

Based on the loss exceedance function, the losses l_t with a specific return period t can be found as

$$\begin{array}{}\text{(2)}& l_{\mathrm{t}}={\mathit{\lambda }}_{\mathrm{L}}^{-\mathrm{1}}\left({\displaystyle \frac{\mathrm{1}}{t}}\right)=F_{\mathrm{L}}^{-\mathrm{1}}\left(\mathrm{1}-{\displaystyle \frac{\mathrm{1}}{{\mathit{\lambda }}_{\mathrm{H}}t}}\right),\end{array}$$

which is defined only for $t\ge \frac{\mathrm{1}}{{\mathit{\lambda }}_{\mathrm{H}}}$. The loss l_t is also called the t-year loss.

Following Rosero-Velásquez and Straub (2022), the representative earthquake scenario θ_t, associated with a return period t, is defined as the most likely scenario among those causing the t-year loss l_t. In other words, θ_t is the mode of the conditional PDF of Θ given the loss L=l_t, also called the loss disaggregation of Θ given L=l_t:

$$\begin{array}{}\text{(3)}& {\mathit{\theta }}_{\mathrm{t}}=\underset{\mathit{\theta }}{\text{arg\hspace{0.17em}\hspace{0.17em}max}}{f}_{\mathbf{\Theta }\mathrm{|}L}\left(\mathit{\theta }\mathrm{|}{l}_{\mathrm{t}}\right).\end{array}$$

wherein $\underset{(\cdot )}{\text{arg\hspace{0.17em}max}(\cdot )}$ returns the point at which the objective function has a maximum. Equation (3) defines the representative earthquake scenario by conditioning on the occurrence of the loss l_t, whereby l_t is defined in terms of the exceedance rate. The equation describes the scenario that is most likely to lead to the t-year loss l_t.

An alternative definition can be formulated in terms of loss exceedance instead of loss occurrence:

$$\begin{array}{}\text{(4)}& {\mathit{\theta }}_{\mathrm{t}}^{\text{exc}}=\underset{\mathit{\theta }}{\text{arg\hspace{0.17em}\hspace{0.17em}max}}{f}_{\mathbf{\Theta }\mathrm{|}L}(\mathit{\theta }\mathrm{|}L\ge l_{\mathrm{t}}).\end{array}$$

Equation (4) defines the scenario that is most likely to exceed l_t. This is the definition corresponding to the classical loss disaggregation (Goda and Hong, 2009; Jayaram and Baker, 2009b). We note that with this definition, in general, the scenario representative of a t-year loss will have a return period higher than t. Hence, we find its interpretation more difficult and prefer the definition in Eq. (3). Nevertheless, we propose algorithms to evaluate the representative scenarios according to the two definitions and compare the resulting scenarios in an illustrative example.

Figure 1 illustrates the conditional distributions f_Θ|L(θ|l_t) and f_Θ|L(θ|L≥l_t) for two source parameters $\mathbf{\Theta }=[{\mathrm{\Theta }}_{\mathrm{1}},{\mathrm{\Theta }}_{\mathrm{2}}]$, e.g., representing the magnitude M_w and the hypocentral distance R with respect to a location of interest.

By Bayes' rule, Eq. (3) can be expressed in terms of f_Θ(θ),

$$\begin{array}{}\text{(5)}& {\mathit{\theta }}_{\mathrm{t}}=\underset{\mathit{\theta }}{\text{arg\hspace{0.17em}\hspace{0.17em}max}}{f}_{L\mathrm{|}\mathbf{\Theta }}\left(l_{\mathrm{t}}\mathrm{|}\mathit{\theta }\right)f_{\mathbf{\Theta }}\left(\mathit{\theta }\right),\end{array}$$

and similarly for the loss exceedance approach,

$$\begin{array}{}\text{(6)}& {\mathit{\theta }}_{\mathrm{t}}^{\text{exc}}=\underset{\mathit{\theta }}{\text{arg\hspace{0.17em}\hspace{0.17em}max}}\left(\mathrm{1}-F_{L\mathrm{|}\mathbf{\Theta }}\left(l_{\mathrm{t}}\mathrm{|}\mathit{\theta }\right)\right)f_{\mathbf{\Theta }}\left(\mathit{\theta }\right),\end{array}$$

wherein f_L|Θ(l|θ) and F_L|Θ(l_t|θ) are, respectively, the conditional PDF and CDF, illustrating that the scenario selection criterion balances the probability of the earthquake scenario (quantified by f_Θ(θ)) and the probability of the t-year losses occurring (or being exceeded) at that scenario.

To ease the notation in the following section, we let z_t(θ) denote the objective function of Eq. (5),

$$\begin{array}{}\text{(7)}& z_{\mathrm{t}}\left(\mathit{\theta }\right)=f_{L\mathrm{|}\mathbf{\Theta }}\left(l_{\mathrm{t}}\mathrm{|}\mathit{\theta }\right)f_{\mathbf{\Theta }}\left(\mathit{\theta }\right),\end{array}$$

and $z_{\mathrm{t}}^{\text{exc}}\left(\mathit{\theta }\right)$ denote the objective function of Eq. (6),

$$\begin{array}{}\text{(8)}& z_{\mathrm{t}}^{\text{exc}}\left(\mathit{\theta }\right)=\left(\mathrm{1}-F_{L\mathrm{|}\mathbf{\Theta }}\left(l_{\mathrm{t}}\mathrm{|}\mathit{\theta }\right)\right)f_{\mathbf{\Theta }}\left(\mathit{\theta }\right).\end{array}$$

3 Method for using scenario selection based on a synthetic earthquake catalog

We consider the case where the randomness of earthquake events is represented through a synthetic earthquake catalog. Specifically, we aim at identifying the earthquake scenarios in the catalog that maximize the objective function in Eqs. (5) and (6) for different return periods t.

The objective functions of Eqs. (5) and (6) consist of the PDF f_Θ(θ), which is known from the earthquake source model, and the conditional PDF or CDF of L given Θ evaluated at l_t, which can be approximated with conditional samples of losses. To account for the aleatory uncertainty in the modeled ground motions one can draw Monte Carlo samples from the catalog (Silva, 2016) and propagate them to the loss metrics. However, performing this number of loss evaluations for an entire seismic catalog (normally containing dozens of thousands of events) is computationally (too) expensive. As an alternative, one can use Gaussian process models in combination with active learning to handle aleatory uncertainty more efficiently (Tomar and Burton, 2021; Rosero-Velásquez and Straub, 2022). Furthermore, pre-selecting scenarios by the use of extreme value theory and the generalized Pareto distribution (Borzoo et al., 2021) has been proposed.

We propose first performing only one loss evaluation for each scenario in the catalog and using these to approximate the loss exceedance function and l_t. The same samples are used for an initial approximation of f_L|Θ(l_t|θ), the second part of the objective function. This approximation is improved by the use of active learning to identify earthquake scenarios in the catalog for which additional loss evaluations are to be performed.

This methodology is an adaptation of the one proposed in Rosero-Velásquez and Straub (2022).

Figure 2 General procedure for selecting representative earthquake scenarios with a synthetic earthquake catalog in terms of loss occurrence. Section 3.1 explains panels (e) and (f) in more detail, and Sects. 3.3 and 3.4 explain panel (h). The remaining panels are referred to in Sect. 3. The procedure in terms of loss exceedance only differs in panel (f), and it is explained in Sect. 3.2.

Figure 2 illustrates the main steps of the methodology for selecting representative earthquake scenarios for n_t return periods $t_{\mathrm{1}}>t_{\mathrm{2}}\mathrm{\dots }>t_{n_{\mathrm{t}}}$. The earthquake model in this example is a single seismic source within a bounding volume, and the system is a single building.

The starting point is a seismic source model (Fig. 2a), which consists of the occurrence rate λ_H and the PDF of the source parameters f_Θ(θ), together with an associated stochastic seismic catalog (Fig. 2b). The catalog consists of a set of n random and independent earthquake scenarios ${\mathit{\theta }}^{\left(\mathrm{1}\right)},\mathrm{\dots },{\mathit{\theta }}^{\left(n\right)}$ generated from f_Θ(θ), possibly associated with weights ${\mathit{\omega }}_{\mathrm{1}},\mathrm{\dots },{\mathit{\omega }}_n$ with ${\sum }_{i=\mathrm{1}}^n{\mathit{\omega }}_i=\mathrm{1}$. The generation of such catalogs for the study site is described in Sect. 5.2.

For each scenario, one simulates the ground motion fields in terms of the intensity measure (e.g., the peak ground acceleration, PGA) through the GMM. These intensity measures are the input to assess the performance of the system components by combining them with vulnerability functions. Based on the component performances, the total losses in the system l are evaluated (Fig. 2c). Details on the simulation of the ground motion, the system response, and loss calculation for the study site are given in Sect. 5.

From these samples of the system losses, one obtains an estimate of the loss exceedance curve ${\widehat{\mathit{\lambda }}}_{\mathrm{L}}\left(l\right)$:

$$\begin{array}{}\text{(9)}& {\widehat{\mathit{\lambda }}}_{\mathrm{L}}\left(l\right)={\mathit{\lambda }}_{\mathrm{H}}\sum _{i=\mathrm{1}}^n{\mathit{\omega }}_i\mathfrak{1}(l^{\left(i\right)}>l),\end{array}$$

where $\mathfrak{1}(\cdot )$ is the indicator function. In addition, one obtains estimates of the t-year losses ${\widehat{l}}_{\mathrm{t}}$ for all return periods of interest $t_{\mathrm{1}},\mathrm{\dots },t_{n_{\mathrm{t}}}$ following Eq. (2) (Fig. 2d).

Since the conditional density of the losses f_L|Θ(l|θ) is not available in analytical form, we propose approximating it with ${\widehat{f}}_{L\mathrm{|}\mathbf{\Theta }}\left(l\mathrm{|}\mathit{\theta }\right)$ (Fig. 2e and f), as detailed in Sect. 3.1. We utilize this approximation in the objective function of Eq. (5) (Fig. 2g) to obtain initial estimates of the objective function z_t(θ) at each scenario of the catalog, which we denote as $z_{\mathrm{t}}^{\left(\mathrm{1}\right)},\mathrm{\dots },z_{\mathrm{t}}^{\left(n\right)}$. To reduce the scatter in the estimates of z_t, we add a smoothing step (Fig. 2h), which is described in Sect. 3.3.

At this and any later stage of the algorithm, we approximate the solution of Eq. (5) by Fig. 2i:

$$\begin{array}{}\text{(10)}& {\displaystyle }& {\displaystyle }{i}_{\mathrm{t}}=\underset{i=\mathrm{1},\mathrm{\dots },n}{\text{arg\hspace{0.17em}\hspace{0.17em}max}}{z}_{\mathrm{t}}^{\left(i\right)},\text{(11)}& {\displaystyle }& {\displaystyle }{\mathit{\theta }}_{\mathrm{t}}\approx {\mathit{\theta }}^{\left(i_{\mathrm{t}}\right)}.\end{array}$$

The initial approximation based on a single loss evaluation per scenario in the catalog is typically poor. To enhance the accuracy, we use an active-learning strategy (Fig. 2h). It intelligently selects earthquake scenarios from the catalog, for which additional loss evaluations are performed. This is presented in Sect. 3.4.

For the representative earthquake scenario defined by the loss exceedance approach, we approximate the conditional CDF of the losses with the empirical CDF ${\widehat{F}}_{L\mathrm{|}\mathbf{\Theta }}\left(l\mathrm{|}\mathit{\theta }\right)$ (analogous to Fig. 2f) and utilize this approximation in the objective function of Eq. (6) (analogous to Fig. 2g) to obtain initial estimates of $z_{\mathrm{t}}^{\text{exc}}\left(\mathit{\theta }\right)$, denoted by $z_{\mathrm{t}}^{\text{exc}\left(\mathrm{1}\right)},\mathrm{\dots },z_{\mathrm{t}}^{\text{exc}\left(n\right)}$. We then approximate the solution of Eq. (6) by

$$\begin{array}{}\text{(12)}& {\displaystyle }& {\displaystyle }{i}_{\mathrm{t}}^{\text{exc}}=\underset{i=\mathrm{1},\mathrm{\dots },n}{\text{arg\hspace{0.17em}\hspace{0.17em}max}}{z}_{\mathrm{t}}^{\text{exc}\left(i\right)},\text{(13)}& {\displaystyle }& {\displaystyle }{\mathit{\theta }}_{\mathrm{t}}^{\text{exc}}\approx {\mathit{\theta }}^{\left(i_{\mathrm{t}}^{\text{exc}}\right)}.\end{array}$$

3.1 Approximation of the objective function z_t(θ) with kernel density estimation

We approximate the conditional density f_L|Θ(l|θ) using weighted kernel density estimation (KDE) (Gisbert, 2003). The KDEs at each scenario θ⁽ⁱ⁾ are evaluated with n₂ loss evaluations, which come from the closest scenarios ${\mathit{\theta }}^{\left(c_{i,\mathrm{1}}\right)},\mathrm{\dots },{\mathit{\theta }}^{\left(c_{i,n_{\mathrm{2}}}\right)}$ and have associated weights $w_{i,\mathrm{1}},\mathrm{\dots },w_{i,n_{\mathrm{2}}}$ which sum up to 1, i.e., ${\sum }_{j=\mathrm{1}}^{n_{\mathrm{2}}}{w}_{i,j}=\mathrm{1}$:

$$\begin{array}{}\text{(14)}& {\widehat{f}}_{L\mathrm{|}\mathbf{\Theta }}\left(l\mathrm{|}{\mathit{\theta }}^{\left(i\right)}\right)=\sum _{j=\mathrm{1}}^{n_{\mathrm{2}}}{w}_{i,j}\mathit{\kappa }\left(l,l^{\left(c_{i,j}\right)},\mathit{\gamma }\right),\end{array}$$

where κ is a kernel function and γ is the bandwidth. We define the weights as $w_{i,j}=\exp (-d_{i,j})/{\sum }_{k=\mathrm{1}}^{n_{\mathrm{2}}}\exp (-d_{i,k})$, where d_i,j is the Mahalanobis distance between θ⁽ⁱ⁾ and ${\mathit{\theta }}^{\left(c_{i,j}\right)}$. This ensures that the loss values from scenarios similar to θ⁽ⁱ⁾ are given more weight in the KDE.

A common choice for κ is the Gaussian kernel function, which employs the standard Gaussian PDF ϕ(⋅),

$$\begin{array}{}\text{(15)}& \mathit{\kappa }\left(l,l^{\left(c_{i,j}\right)},\mathit{\gamma }\right)={\displaystyle \frac{\mathrm{1}}{\mathit{\gamma }}}\mathit{\varphi }\left({\displaystyle \frac{l-l^{\left(c_{i,j}\right)}}{\mathit{\gamma }}}\right),\end{array}$$

with γ computed as suggested in Silverman (1986). Alternatively, one can employ a lognormal kernel, excluding the zero loss values (if any), whose probability $p_{\mathrm{0}}^{\left(i\right)}$ is estimated from the conditional-loss samples $l^{\left(c_{i,\mathrm{1}}\right)},\mathrm{\dots },l^{\left(c_{i,n_{\mathrm{2}}}\right)}$. That is, for l>0 and $l^{\left(c_{i,j}\right)}>\mathrm{0}$,

$$\begin{array}{}\text{(16)}& \mathit{\kappa }\left(l,l^{\left(c_{i,j}\right)},\mathit{\gamma }\right)=\left(\mathrm{1}-p_{\mathrm{0}}^{\left(i\right)}\right){\displaystyle \frac{\mathrm{1}}{l\mathit{\gamma }}}\mathit{\varphi }\left({\displaystyle \frac{\ln l-\ln l^{\left(c_{i,j}\right)}}{\mathit{\gamma }}}\right),\end{array}$$

wherein the bandwidth γ is computed as suggested in Silverman (1986) but only using the logarithm of the nonzero loss samples. As a consequence, the weights w_i,j have to be adjusted excluding the zero loss samples, i.e.,

$$\begin{array}{}\text{(17)}& w_{i,j}=\mathfrak{1}\left(l^{\left(c_{i,j}\right)}>\mathrm{0}\right){\displaystyle \frac{\exp \left(-d_{i,j}\right)}{{\sum }_{k=\mathrm{1}}^{n_{\mathrm{2}}}\mathfrak{1}\left(l^{\left(c_{i,k}\right)}>\mathrm{0}\right)\exp \left(-d_{i,k}\right)}}.\end{array}$$

In this case, for scenarios where all the n₂ conditional-loss samples are zero, the density at l_t equals zero.

The choice of n₂ for the KDEs is associated with a trade-off: on the one hand, a small n₂ leads to a poor density estimation but one that is based on loss samples coming from similar scenarios. On the other hand, a large n₂ produces a biased KDE from the true conditional density, since it incorporates loss samples of more dissimilar scenarios. However, the bias can be reduced by additional model evaluations. In fact, n₂ or more model evaluations at a scenario θ provide a more accurate KDE than a KDE based on model evaluations coming from the n₂ closest scenarios to θ.

For each return period, we obtain an estimate of f_L|Θ(l_t|θ⁽ⁱ⁾) by evaluating Eq. (14) with the argument ${\widehat{l}}_{\mathrm{t}}$. This process is illustrated in Fig. 2e and f. By multiplication with the prior, an estimate of the objective function at all scenarios in the catalog is obtained:

$$\begin{array}{}\text{(18)}& z_{\mathrm{t}}^{\left(i\right)}={\widehat{f}}_{L\mathrm{|}\mathbf{\Theta }}\left({\widehat{l}}_{\mathrm{t}}\mathrm{|}{\mathit{\theta }}^{\left(i\right)}\right)f_{\mathbf{\Theta }}\left({\mathit{\theta }}^{\left(i\right)}\right).\end{array}$$

To reduce the significant noise associated with these estimates, we update them with an additional smoothing step described in Sect. 3.3.

Even after the additional smoothing step, the estimates $z_{\mathrm{t}}^{\left(i\right)}$ remain subject to uncertainty, due to the limited number of noisy loss evaluations and the need to pool the evaluations from multiple scenarios. For the purpose of the active-learning procedure presented in Sect. 3.4, we approximate the uncertainty associated with the objective function values by modeling the estimates $z_{\mathrm{t}}^{\left(i\right)}$ as Gaussian random variables $Z_{\mathrm{t}}^{\left(i\right)}$. We denote their mean values as ${\mathit{\mu }}_{Z_{\mathrm{t}}}^{\left(i\right)}$ and their standard deviations as ${\mathit{\sigma }}_{Z_{\mathrm{t}}}^{\left(i\right)}$. The mean values are set to $z_{\mathrm{t}}^{\left(i\right)}$. We estimate the standard deviations ${\mathit{\sigma }}_{Z_{\mathrm{t}}}^{\left(i\right)}$ for $i=\mathrm{1},\mathrm{\dots },n$ via bootstrapping (Efron and Tibshirani, 1993) on their n₂ nearest-neighbors n_b times, wherein conditional losses are resampled according to the weights w_i,j.

3.2 Approximation of the objective function $z_{\mathrm{t}}^{\text{exc}}\left(\mathit{\theta }\right)$ with a weighted empirical CDF

Equation (6) contains the conditional CDF of the losses given a scenario. It can be approximated at θ⁽ⁱ⁾ with a weighted empirical CDF based on the n₂ loss evaluations coming from the closest scenarios and their associated weights:

$$\begin{array}{}\text{(19)}& {\widehat{F}}_{L\mathrm{|}\mathbf{\Theta }}\left(l\mathrm{|}{\mathit{\theta }}^{\left(i\right)}\right)=\sum _{j=\mathrm{1}}^{n_{\mathrm{2}}}{w}_{i,j}\mathfrak{1}(l^{\left(c_{i,j}\right)}<l).\end{array}$$

The objective function evaluated at scenario i is then estimated as follows:

$$\begin{array}{}\text{(20)}& z_{\mathrm{t}}^{\text{exc}\left(i\right)}=\left(\mathrm{1}-{\widehat{F}}_{L\mathrm{|}\mathbf{\Theta }}\left({\widehat{l}}_{\mathrm{t}}\mathrm{|}{\mathit{\theta }}^{\left(i\right)}\right)\right)f_{\mathbf{\Theta }}\left({\mathit{\theta }}^{\left(i\right)}\right).\end{array}$$

We apply to these estimates the same uncertainty treatment described in Sect. 3.1 for the KDEs. Thus, we model the estimates as Gaussian random variables $Z_{\mathrm{t}}^{\text{exc}\left(i\right)}$ with a mean ${\mathit{\mu }}_{Z_{\mathrm{t}}^{\text{exc}}}^{\left(i\right)}=z_{\mathrm{t}}^{\text{exc}\left(i\right)}$ and standard deviation ${\mathit{\sigma }}_{Z_{\mathrm{t}}^{\text{exc}}}^{\left(i\right)}$ estimated via bootstrapping.

3.3 Smoothed estimation of the objective function with Gaussian process regression

To reduce the noise in the estimates of the objective function, we perform an additional smoothing step via Gaussian process regression (GPR) (Rasmussen and Williams, 2005). For each return period t and in each step of the active-learning algorithm described in Sect. 3.4, we perform a separate GPR.

A drawback of GPR is that the computational cost escalates with the size of the training set n_train. Fitting and estimating the objective function using standard GPR is an 𝒪(n³) task (Rasmussen and Williams, 2005). Therefore, we perform GPR smoothing only for estimates θ⁽ⁱ⁾ near the current solution of Eq. (12), and the GPR hyperparameters are learned only once in the first step. Specifically, only for the training set do we identify the n_train=1500 nearest scenarios using the Mahalanobis distance, train the GPR, and replace the estimate of $z_{\mathrm{t}}^{\left(i\right)}$ ($z_{\mathrm{t}}^{\text{exc}\left(i\right)}$). The other estimates of $z_{\mathrm{t}}^{\left(i\right)}$ ($z_{\mathrm{t}}^{\text{exc}\left(i\right)}$) are left unaltered.

3.4 Active learning

An accurate estimation of the objective function is only important near the solution. We exploit this by employing an active-learning (AL) strategy to identify scenarios for which further model evaluations are performed.

AL selects scenarios to evaluate through the acquisition function. Here we use the augmented expected improvement (AEI) as an acquisition function (Huang et al., 2006). It approximates for each scenario the expected value of the improvement of the objective function over the current maximum.

We modify the AEI of Huang et al. (2006) with a correction factor, which assesses the quality of the KDE at each scenario. The resulting AEI at scenario θ⁽ⁱ⁾ is

$$\begin{array}{}\text{(21)}& \mathrm{AEI}\left({\mathit{\theta }}^{\left(i\right)}\right)={\displaystyle \frac{\mathrm{1}}{c_{\text{neigh}}^{\left(i\right)}}}\mathbb{E}\left[max\left(Z_{\mathrm{t}}^{\left(i\right)}-z_{\mathrm{t}}^{\ast },\mathrm{0}\right)\right],\end{array}$$

wherein $z_{\mathrm{t}}^{\ast }$ is the estimate of the objective function at the current best solution θ^∗, which is defined as (Huang et al., 2006)

$$\begin{array}{}\text{(22)}& {\mathit{\theta }}^{\ast }=\underset{i=\mathrm{1},\mathrm{\dots },n}{\text{arg\hspace{0.17em}\hspace{0.17em}max}}\left(z_{\mathrm{t}}^{\left(i\right)}+{\mathit{\sigma }}_{Z_{\mathrm{t}}}^{\left(i\right)}\right).\end{array}$$

The factor $c_{\text{neigh}}^{\left(i\right)}$ considers the KDE estimation quality at θ⁽ⁱ⁾. We define it as follows:

$$\begin{array}{}\text{(23)}& c_{\text{neigh}}^{\left(i\right)}=max\left(\sum _{j=\mathrm{1}}^{n_{\mathrm{2}}}\exp \left(-d_{i,j}\right),n^{\left(i\right)}\right),\end{array}$$

wherein n⁽ⁱ⁾ is the sample size of conditional-loss values simulated at θ⁽ⁱ⁾.

The expected value in Eq. (21) is computed in terms of the standard normal PDF ϕ(⋅) and CDF Φ(⋅) (Huang et al., 2006):

$$\begin{array}{}\text{(24)}& \begin{array}{rl}\mathbb{E}\left[max\left(Z_{\mathrm{t}}^{\left(i\right)}-z_{\mathrm{t}}^{\ast },\mathrm{0}\right)\right]& =\left({\mathit{\mu }}_{Z_{\mathrm{t}}}^{\left(i\right)}-z_{\mathrm{t}}^{\ast }\right)\mathrm{\Phi }\left({\displaystyle \frac{{\mathit{\mu }}_{Z_{\mathrm{t}}}^{\left(i\right)}-z_{\mathrm{t}}^{\ast }}{{\mathit{\sigma }}_{Z_{\mathrm{t}}}^{\left(i\right)}}}\right)\\ & +{\mathit{\sigma }}_{Z_{\mathrm{t}}}^{\left(i\right)}\mathit{\varphi }\left({\displaystyle \frac{{\mathit{\mu }}_{Z_{\mathrm{t}}}^{\left(i\right)}-z_{\mathrm{t}}^{\ast }}{{\mathit{\sigma }}_{Z_{\mathrm{t}}}^{\left(i\right)}}}\right).\end{array}\end{array}$$

For each return period t, we perform n_l loss evaluations at the n_s scenarios with the largest AEI. Taking into account the $n_{\mathrm{s}}\times n_{\mathrm{l}}\times n_{\mathrm{t}}$ new model evaluations, we update the KDEs, the density observations $z_{\mathrm{t}}^{\left(i\right)}$, and the bootstrap standard deviations ${\mathit{\sigma }}_{Z_{\mathrm{t}}}^{\left(i\right)}$. At scenarios where more than n₂ loss evaluations have been computed, we deviate from Eq. (14) and evaluate the KDE with all these evaluations (instead of only n₂ evaluations).

The AL steps are repeated until convergence is achieved or the maximum number of AL iterations n₃ is exceeded. Convergence is achieved when the AEI of all scenarios is below a threshold ϵ for at least n_d consecutive AL iterations, which prevents premature stopping. A suggested value for n_d is d+1, wherein d is the dimensionality of the source parameter random vector Θ (Huang et al., 2006). We also choose n₃=1000 for encouraging the AL procedure to stop by convergence. The threshold ϵ is chosen as (Huang et al., 2006)

$$\begin{array}{}\text{(25)}& \mathit{ϵ}=r\times \left(\underset{i=\mathrm{1},\mathrm{\dots },n}{max}{z}_{\mathrm{t}}^{\mathrm{0},\left(i\right)}-\underset{i=\mathrm{1},\mathrm{\dots },n}{min}{z}_{\mathrm{t}}^{\mathrm{0},\left(i\right)}\right),\end{array}$$

where $z_{\mathrm{t}}^{\mathrm{0},\left(i\right)}$ is the initial evaluation of the objective function, which is computed before the AL.

An analogous derivation of the AEI is obtained for the case of the objective function in terms of loss exceedance, i.e., $z_{\mathrm{t}}^{\text{exc}\left(i\right)}$.

4 Illustrative examples

In this section, we present two simple examples to illustrate the methodology. The first one is a one-dimensional example, where we focus on the performance of AL and the approximations of the objective function obtained with the noisy KDE estimations and GPR. The second one is a two-dimensional example, where we show the variability of the scenario selection and the solutions computed with the loss occurrence and exceedance approaches. In both examples, the exact solution is known.

Figure 3 Illustration of active learning (AL) for maximizing the objective function z_t(m_w) with the fixed value l_t, marked with a dashed line in the plots in the left panel (a, c). The solution m_w,t is approximated with the sample point ${\widehat{m}}_{w,t}$. The loss samples before and during the AL steps are shown in the left panel (b, d). The approximations of the objective function, either with kernel density estimation (KDE) or Gaussian process regression (GPR), are shown in the right panel (b, d).

4.1 Residential-building stock subjected to a single seismic source with a variable magnitude

Figure 3 illustrates the performance of the acquisition function on a one-dimensional example. In the example, the only source parameter is the magnitude M_w, which follows beta distribution with shape parameters α=1 and β=3 that is scaled between 0 and 10. The conditional distribution of the logarithm of the losses associated with the residential-building stock given M_w=m_w is a normal random variable with a mean ${\mathit{\mu }}_{\ln L}\left(m_{\mathrm{w}}\right)=-\frac{\mathrm{1}}{\mathrm{2}}\sin \left(\frac{\mathrm{5}}{\mathrm{2}}\right(m_{\mathrm{w}}-\mathrm{5}\left)\right)+\exp \left(\frac{m_{\mathrm{w}}-\mathrm{5}}{\mathrm{7}}\right)+\mathrm{7}$ and standard deviation σ_ln L=0.7. We set l_t=10 for this example.

We generate a synthetic earthquake catalog of size n=100. The KDEs are computed with the Gaussian kernel and n₂=70 loss samples or more if the scenario has more than n₂ conditional-loss evaluations, and the bootstrap standard deviation is computed based on n_b=100 samples. We perform the GPR on the whole catalog and learn the hyperparameters at every AL step, since the catalog size in this example is not restrictive. For the AL stage, we select n_s=5 scenarios per AL iteration for computing n_l=10 loss evaluations at each scenario (i.e., $n_{\mathrm{s}}\times n_{\mathrm{l}}=\mathrm{50}$ damage evaluations per AL iteration). The acquisition function is the AEI, as introduced in Eq. (21). We also let the algorithm achieve convergence with a maximum of n₃=1000 AL iterations, with the convergence criterion in Eq. (25) and r=0.001.

Figure 3 compares intermediate and final results for this example to the true results. After the initial loss evaluations at the 100 scenarios, the estimate of the objective function is poor. However, the acquisition function is able to select scenarios near the true solution. In the final step, one can observe that estimates of the objective function values have high noise, but the GPR is effective in reducing this noise. The resulting estimate of the objective function is close to the true value around the optimum.

4.2 Residential-building stock subjected to an earthquake with an unknown magnitude and location

This simple example is adapted from Rosero-Velásquez and Straub (2022). It considers a hypothetical fault, where strong earthquakes occur with a rate of ${\mathit{\lambda }}_{\mathrm{H}}=\mathrm{0.3}{\text{yr}}^{-\mathrm{1}}$. We consider the damage that earthquakes cause to a residential-building stock in a small town. The source parameters $\mathbf{\Theta }=[M_{\mathrm{w}},\ln R{]}^{⊺}$ are the magnitude M_w and the average hypocentral distance R from the earthquake source to the buildings. The source model for Θ (f_Θ(θ)) is a normal distribution with a mean vector μ_Θ and covariance matrix Σ_Θ given as follows:

$$\begin{array}{}\text{(26)}& {\mathit{\mu }}_{\mathbf{\Theta }}=\left[\begin{array}{c}\mathrm{7.00}\\ \mathrm{4.38}\end{array}\right],{\mathbf{\Sigma }}_{\mathbf{\Theta }}=\left[\begin{array}{cc}\mathrm{0.36}& -\mathrm{0.08}\\ -\mathrm{0.08}& \mathrm{0.49}\end{array}\right].\end{array}$$

A standard deviation σ represents the uncertainty in the ground motion, damage measure, and losses. The losses L are a log-normal random variable with parameters ${\mathit{\mu }}_{\ln L}=-\mathrm{3.16}$ and ${\mathit{\sigma }}_{\ln L}=\sqrt{\mathrm{2.46}+{\mathit{\sigma }}^{\mathrm{2}}}$. With these choices, the conditional density f_Θ|L(θ|l_t) can be evaluated analytically. It is a normal distribution, whose mean vector is the representative earthquake scenario for a return period t.

We set σ=0.5 and use return periods of 50, 100, 500, and 1000 years. The resulting exact representative earthquake scenarios with the loss occurrence approach are ${\mathit{\theta }}_{\mathrm{50}}=[\mathrm{7.42},\mathrm{3.42}{]}^{⊺}$, ${\mathit{\theta }}_{\mathrm{100}}=[\mathrm{7.51},\mathrm{3.21}{]}^{⊺}$, ${\mathit{\theta }}_{\mathrm{500}}=[\mathrm{7.69},\mathrm{2.81}{]}^{⊺}$, and ${\mathit{\theta }}_{\mathrm{1000}}=[\mathrm{7.75},\mathrm{2.65}{]}^{⊺}$. We use them to verify the proposed sampling-based algorithm.

To estimate θ_t with the proposed methodology, a synthetic earthquake catalog with $n=\mathrm{2}\times {\mathrm{10}}^{\mathrm{4}}$ random scenarios is employed. We simulate the losses once at each scenario, approximate l_t, and compute the KDEs at each scenario. The KDEs are based on n₂=200 loss values, computed with the Gaussian kernel, and the bootstrap variance is computed with n_b=100 repetitions. We perform GPR with a training set of size n_train=1500, which is constructed as described in Sect. 3.3.

Figure 4 Numerical approximation of the representative earthquake scenarios. Panel (a) shows the representative scenarios computed with the loss occurrence approach ${\widehat{\mathit{\theta }}}_{\mathrm{t}}$, and panel (b) shows those computed with the loss exceedance approach ${\widehat{\mathit{\theta }}}_{\mathrm{t}}^{\text{exc}}$. The representative earthquake scenarios correspond to four different return periods, t=50, 100, 500, and 1000 years, based on a Monte Carlo sample of scenarios. Each return period is represented by a different color. For each return period, the 20 approximations of ${\widehat{\mathit{\theta }}}_{\mathrm{t}}$ (${\widehat{\mathit{\theta }}}_{\mathrm{t}}^{\text{exc}}$), corresponding to 20 experiments, are the colored empty circles, and the corresponding exact solutions are depicted by filled circles. The grey points are the scenarios of the catalog, and the dashed contours represent the probability density function (PDF) of the source parameters.

The maximum number of AL iterations is n₃=1000, where at every step the losses are evaluated at n_s=2 scenarios n_l=10 times, for each return period. The procedure stops after the maximum AEI is below ϵ, with r=0.001, for at least n_d=5 consecutive AL iterations. For analyzing the uncertainty in the estimation of θ_t, we repeat the experiment 20 times. For the 50-, 100-, and 500-year representative scenarios, all experiments converged in less than 10 AL iterations, whereas for the 1000-year return period at most 30 AL iterations were required. For the loss exceedance approach, fewer iterations were required in general.

Figure 4 shows the resulting representative hazard scenarios for each return period and their spread, which is mainly caused by the numerical approximation of the objective function with a limited number of samples. As expected, one can observe that the representative scenarios are more extreme when using the loss exceedance approach.

5 Case study: the communes of Valparaíso and Viña del Mar

5.1 Context of the study area

We apply the proposed methodology to determine representative earthquake scenarios for the communes of Valparaíso and Viña del Mar, which are located in the Chilean region of Valparaíso on the Pacific coast. The study area is the second-largest Chilean urban center; based on the latest Chilean census (INE, 2017), it is home to 630 903 inhabitants. It hosts the Port of Valparaíso, which is an important container port and the main passenger port in Chile. The area has a heterogeneous building inventory, ranging from apartment buildings to informal settlements and a historic district declared a World Heritage Site by UNESCO in 2003 (Indirli et al., 2011; Jiménez et al., 2018).

The National Electric System (SEN) provides the area's power supply. The SEN is the largest Chilean transmission grid and covers most of the national territory. The SEN connects power plants and substations with the consumer areas through high-voltage lines. The topology of the SEN is characterized as a single-scale network with a fast decaying tail, and most of the load substations are close to a generation unit, with a median distance of 9 km (Ferrario et al., 2022).

Powerful earthquakes have hit the area in the past, such as the 1730 earthquake, with an inferred magnitude M_w in the range of 9.1–9.3, and the 1906 event, with an inferred moment magnitude M_w of 8.0–8.2 (Carvajal et al., 2017). More recently, the 1985 M_w 8.0 event affected around 230 000 dwellings and 1 million people in the regions of Valparaíso, Maule, O'Higgins, and the Santiago Metropolitan Region and caused residential-building losses of about USD 1.4 billion (in 1985 values) (ONEMI, 1985). The most recent M_w 8.8 Maule earthquake (2010) caused severe structural damage in buildings in Viña del Mar, including in buildings retrofitted in 1985 (Jünemann et al., 2015).

Table 1 Representative earthquake scenarios for the urban area of Valparaíso selected by Indirli et al. (2011) from historic events. The epicenter location is reported with a map by Indirli et al. (2011), hence their numeric values; the moment magnitude and depth are reproduced here from the earthquake records of the USGS ComCat (Comprehensive Earthquake Catalog) (USGS, 2023). The 1985 M_w 8.0 earthquake event (33°14^′ S, 71°51^′ W; depth of 33 km) has source parameters similar to those of the 1906 event, while the 2010 earthquake event occurred after the study of Indirli et al. (2011) was submitted for publication. The event dates are in local time. The locations of the epicenters are shown in the map of Fig. 11.

A hazard evaluation of the Valparaíso urban area presented by Indirli et al. (2011) selected representative earthquake scenarios based on historical events, considering the seismicity around the study area, to specify an average regional seismic input and to generate synthetic seismograms. The scenarios are summarized in Table 1; their magnitudes range from M_w=5.7 to M_w=8.2.

5.2 Earthquake model and synthetic earthquake catalog

We employ the earthquake model presented by Poulos et al. (2019) to generate a synthetic catalog of earthquake scenarios. The catalog has 2 × 10⁴ scenarios with a magnitude larger than or equal to M_w=5.0, which is the minimum magnitude defined by Poulos et al. (2019) for performing the declustering on the historical seismic catalogs on which the earthquake model is based on. The catalog covers the whole country of Chile and consists of scenarios at the subduction interface and subduction intraslab zones. The earthquake model utilizes the slab geometry proposed by Hayes et al. (2012) for the depth contours and trench geometry and divides the Chilean subduction zone into three subduction interface and four intraslab zones, whose combined occurrence rate equals λ_H = 43.3 yr⁻¹.

The epicentral locations of the catalog are generated randomly, based on the occurrence rate associated with the seismic zones defined by the occurrence model. The magnitude is sampled with an importance sampling (IS) approach. We employ a uniform distribution with minimum and maximum values defined by the magnitude range of each seismic zone as the IS density. The corresponding IS weights are considered when determining the loss exceedance function and computed following the original Gutenberg–Richter relationship at each seismic zone (Poulos et al., 2019). The resulting catalog is depicted in Fig. 5, in which one can observe that events of different magnitudes have a similar spread within the seven seismic zones.

Figure 5 Synthetic earthquake catalog with 20 000 scenarios (Poulos et al., 2019). The circle size corresponds to the scenario magnitude. The red square contains the study area. Seismic zones 1 to 3 are of the subduction interface type, and zones 4 to 7 are of the subduction intraslab type. Basemap from © OpenStreetMap contributors 2023. Distributed under the Open Data Commons Open Database License (ODbL) v1.0.

The earthquake model of Poulos et al. (2019) only considers the subduction zone, hence the independent source parameters are the moment magnitude M_w, longitude X, and latitude Y of the epicenter. Other parameters, such as the depth H and the strike, dip, and rake angles, are determined by the geometry derived by Hayes et al. (2012), depending of the epicenter location. Therefore, the PDF of the source parameters f_Θ(θ) is represented, respectively, by the conditional PDF $f_{M_{\mathrm{w}}\mathrm{|}X,Y}(m_{\mathrm{w}}\mathrm{|}x,y)$ and the location-dependent occurrence rate λ(x,y):

$$\begin{array}{}\text{(27)}& f_{\mathbf{\Theta }}\left(\mathit{\theta }\right)\propto \mathit{\lambda }(x,y)f_{M_{\mathrm{w}}\mathrm{|}X,Y}(m_{\mathrm{w}}\mathrm{|}x,y).\end{array}$$

Figure 6 (a) Exposure model of the residential-building stock. Each rectangular cell shows the total count of residential buildings, indicating the most dense areas. Source: Pittore et al. (2021). (b) Geographic location of the National Electric System (SEN). The network follows the narrow shape of the country, and the communes of Valparaíso and Viña del Mar (inside red square) are in a central location within the network. Source: Coordinador Eléctrico Nacional (2024), version 2019. Basemap from © OpenStreetMap contributors 2023. Distributed under the Open Data Commons Open Database License (ODbL) v1.0.

5.2.1 Ground motion models

For the residential-building stock, we evaluate the PGA and spectral accelerations at 0.3 and 1.0 s with the ground motion model (GMM) presented by Montalva et al. (2017). The uncertainty in the median prediction is modeled with a Gaussian random field, with the spatial correlation model of Jayaram and Baker (2009a). The choice of these models is based on the epistemic uncertainty analysis of different ground motion and correlation models by Gómez-Zapata et al. (2022a), who analyzed the same study area. For the scope of this study, we do not consider cross-correlated ground motion fields.

For the Chilean power network, we employ the GMM of Abrahamson et al. (2016) and the spatial correlation model developed by Goda and Atkinson (2010) for predicting the PGA. This is the same ground motion model as the one utilized by Ferrario et al. (2022).

The functional form of both GMMs is similar, and therefore, their predictions do not differ significantly, as observed in previous studies (e.g.,  Hussain et al., 2020; Gómez-Zapata et al., 2022a). In particular, Hussain et al. (2020) found negligible differences in direct loss estimates for the residential-building stock of Santiago de Chile after using these two GMMs to simulate the associated ground motion from subduction earthquake scenarios.

5.3 Model for the residential-building stock in the communes of Valparaíso and Viña del Mar

We employ the Bayesian exposure model of the building stock with the building classes described in Gómez-Zapata et al. (2022a) and available in Pittore et al. (2021). The model was constructed by taking the OpenStreetMap footprint of the buildings in the two communes and assigning to each footprint the most likely building class. The buildings are counted within a regular 500 m × 500 m resolution grid in the urban areas, as shown in Fig. 6. Detailed building counts for each class are presented in Gómez-Zapata et al. (2022a)

The model considers 16 building classes (Gómez-Zapata et al., 2022a) which correspond to the ones proposed in the SARA project (Yepes-Estrada et al., 2017) and have an associated replacement cost. Furthermore, each building class has an associated fragility model with five damage states (Villar-Vega et al., 2017). The fragility model for each building associates an intensity measure (spectral acceleration at 0.3 and 1.0 s or the PGA) with the probabilities of achieving a damage state. We assume the following relative replacement cost percentages for each damage level: 0 % for no damage, 2 % for slight damage, 10 % for moderate damage, 50 % for extensive, and 100 % for complete damage.

We utilize the model to evaluate the ground motion and simulate the building damage. Given an intensity level of the ground motion, the damage is simulated randomly at each building with a discrete distribution with probabilities defined by the fragility functions. The losses for each scenario in the catalog are computed as the accumulated reconstruction cost of the damaged residential buildings, based on the simulated damage.

5.4 Model for the Chilean National Electric System (SEN)

We model the SEN and its components following Ferrario et al. (2022). The network model consists of 1494 nodes, representing 500 generation units and 994 substations, as well as the transmission lines connecting them, with a total power generation capacity of 21.9 GW. The model considers seismic interaction and system performance subjected to component failures. Given a scenario with a ground motion field (in this case, the PGA), each node is randomly associated with a damage state, by means of the fragility function, and a recovery time associated with their damage state.

The losses associated with the SEN due to an earthquake scenario are quantified in terms of the energy not supplied (ENS). The ENS is evaluated at each substation by solving the power in normal steady-state operation through the direct current optimal power flow (DCOPF) model (Wood et al., 2013) and comparing it with the power in a damaged state operation caused by the earthquake scenario. To quantify the loss in the power supply in the communes of Valparaíso and Viña del Mar, we calculate the total ENS with the sum of the ENS of all substations located in the two communes (14 in total).

DCOPF is typically adopted in practice for transmission networks (Frank and Rebennack, 2016). It optimizes the power generation cost, taking into account the capacity of the power plants and transmission lines connected to the power grid, the generation cost associated with each power plant, and the demand from the clients. For modeling the system response to an earthquake, the DCOPF model considers the reduced capacity of components affected by the earthquake. A detailed description and validation of the network model of the SEN can be found in Ferrario et al. (2022).

6 Evaluation of representative earthquake scenarios for the communes of Valparaíso and Viña del Mar

6.1 Results for the residential-building stock

Figure 7 shows the annual exceedance rate of the losses (in 2016 USD), based on the replacement costs estimated by Yepes-Estrada et al. (2017). The USGS ComCat records that between 1960 and 2020 there were 12 seismic events that produced a macroseismic intensity greater or equal than VI on the Mercalli scale in the two communes (USGS, 2023). This corresponds to an occurrence rate of 0.20 yr⁻¹. It is reasonable to assume that events with a macroseismic intensity of VI or higher lead to losses of at least USD 10 million. For comparison, the reparation cost of the residential buildings in the two communes due to the 1985 earthquake event was USD 49.7 million (in 1985, which corresponds to USD 110.9 million in 2016) (ODEPLAN, 1985). These data therefore validate the lower end of the loss exceedance rate obtained with the synthetic earthquake catalog.

Figure 7 Loss exceedance function of the reparation costs (USD in 2016 values) associated with the residential-building stock in the communes of Valparaíso y Viña del Mar.

As in Sect. 4.2, we evaluate the representative scenarios in 20 independent runs of the algorithm to check the robustness of the results. In all evaluations, we found a spread of the identified representative scenarios similar to that of Fig. 4. This spread is larger for higher return periods, but most of the numerical solutions (11 out of 20 for the 1000-year loss return period and at least 16 out of 20 for the other loss return periods) have an epicentral location within a radius of 50 km around the mode, and the coefficient of variation in the magnitude is below 4 % for all return periods. In the following, we only present the modes, i.e., the representative scenarios that were identified the most frequently in the 20 repetitions.

Figure 8 shows the representative earthquake scenarios for the analyzed return periods with a loss occurrence approach. One can observe that large return periods are associated with scenarios that have a larger magnitude. The fact that the magnitude for the 1000-year scenario equals only M_w=7.01 is a consequence of the size of the study area. On the one hand, among the 36 seismic events with M_w≥8.0 registered along the Chilean coast between 1570 and 2023 (CSN, 2023), only 4 events had an epicenter near the two communes, i.e., within a radius of approximately 70 km. For comparison, the identified 1000-year scenario is at a distance of around 20 km from Valparaíso and Viña del Mar. On the other hand, the spatial correlation of the ground motion within a small study area leads to an increased likelihood of extreme losses in a scenario with a lower earthquake magnitude and extreme ground motion residuals. This tendency was also found by Goda and Hong (2009) with the loss exceedance approach.

Figure 8 Representative earthquake scenarios for the residential-building stock in the communes of Valparaíso and Viña del Mar for four different return periods (RPs). The hypocentral depth H is displayed in kilometers. Source of the exposure model: Pittore et al. (2021). Basemap from © OpenStreetMap contributors 2023. Distributed under the Open Data Commons Open Database License (ODbL) v1.0.

6.2 Results for the power network

Figure 9 shows the loss exceedance function in terms of the ENS obtained with the synthetic earthquake catalog. The largest sampled ENS value is around 2 × 10⁵ MW h, which is around 20 % of the annual energy demand of the two communes.

Figure 9 Loss exceedance function of the total energy in megawatt hours (MW h) not supplied in communes of Valparaíso and Viña del Mar.

The spread in epicentral locations of the representative scenarios obtained with the 20 runs is larger than the one of the residential-building stock but is still small. At least 13 solutions cluster around the sample mode within a radius of 100 km, and the coefficient of variation in the magnitude is below 5 % for all return periods.

Figure 10 shows the resulting representative earthquake scenarios for the analyzed return periods. One can observe that the scenarios are close to the two communes but less concentrated than those of the residential-building stock and have a different magnitude range. This reflects the fact that the total ENS, although computed only at the substations located within the two communes, depends on the damage state of the components of the rest of network.

Figure 10 Representative earthquake scenarios for the power supply considering the total energy not supplied (ENS) of the communes of Valparaíso and Viña del Mar for four different return periods (RPs). The red square indicates the location of the two communes, and the hypocentral depth H is displayed in kilometers. Source of the power network: Coordinador Eléctrico Nacional (2024). Basemap from © OpenStreetMap contributors 2023. Distributed under the Open Data Commons Open Database License (ODbL) v1.0.

Even though the power supply network is spread out over a larger area, the representative earthquake scenarios are close to the study area. They reflect that the most important components of the network for the two communes are in their proximity. For example, the source location of the 100-year return period scenario lies near a main connection between the substations in the two communes and the rest of the SEN.

6.3 Comparison with past earthquake events

Figure 11 compares the results with the historical events selected in Indirli et al. (2011). Although the representative earthquake scenarios and the selected historical events target the same area of interest, they have different purposes. The historical events presented by Indirli et al. (2011) aim at representing the seismicity of the most important seismic zones affecting the study area. In contrast, the representative earthquake scenarios, as defined in Rosero-Velásquez and Straub (2022), take into account the performance and the losses caused by damage and failures in the analyzed engineering system. In addition, the scenarios are selected based on different loss levels, which are attached to return periods.

Figure 11 Representative earthquake scenarios for the (a) building stock and (b) power supply, compared with past earthquake events selected in Indirli et al. (2011). The location of the two communes is within the red square, and the National Electric System (SEN) network is also displayed. Source of the exposure model: Pittore et al. (2021). Source of the power network: Coordinador Eléctrico Nacional (2024). Basemap from © OpenStreetMap contributors 2023. Distributed under the Open Data Commons Open Database License (ODbL) v1.0.

6.4 Computational cost

In terms of loss evaluations, the analysis requires one evaluation per scenario in the catalog for constructing the loss exceedance function with event-based earthquake risk assessment. That corresponds to 2 × 10⁴ loss evaluations. In addition, during the AL stage, around 10 iterations were necessary to achieve the convergence criterion of Eq. (25), each of them consisting of 160 new loss evaluations (n_s=2 scenarios evaluated n_l=20 times, for each of the n_t=4 return periods). Therefore, 1600 loss evaluations are needed to find the representative earthquake scenarios for four different return periods.

For comparison, Goda and Hong (2009) report that they use a total of 5 × 10⁶ loss evaluations for the classical loss disaggregation. Furthermore, they only evaluate the scenarios with the loss exceedance approach. Extending the loss disaggregation approach to loss occurrence will likely require additional evaluations. Additionally, the computation cost of the loss disaggregation approach scales exponentially with the number of parameters describing seismic scenarios. Hence, the classical loss disaggregation approach will not be applicable to problems in which earthquake scenarios are described by more than three or four parameters. By contrast, we successfully tested the proposed approach for seismic hazard models with seven parameters in applications not reported in this paper.

7 Discussion

7.1 On the results for the communes of Valparaíso and Viña del Mar

The representative earthquake scenarios summarized in Fig. 11 can provide important input to risk assessment and risk management activities. The fact that the scenarios identified with the proposed approach differ from the historical events selected in Indirli et al. (2011) should not be surprising, as the latter are in some sense just “random samples” of earthquake events. Nevertheless, the historic events can provide a useful validation of the identified scenarios. In this regard, the scenarios identified as representative of the power supply network appear to be in line with the historic events. The identified 500- and 1000-year scenarios have larger magnitudes than the historical events, which is expected since the historical events come from a (roughly) 100-year period, as shown in Table 1. By contrast, the representative scenarios identified for the building stock have smaller magnitudes than the historical events. However, they occur much closer to the considered building stock. Furthermore, according to the employed model, extreme losses are more likely to occur through a combination of a less strong earthquake with larger-than-average ground motions (i.e., a large value of the interevent term in the GMM). This effect occurs for the residential-building stock due to its spatial concentration and not (or to a much smaller extent) for the power supply network, which is spatially distributed.

The above observations lead to some important conclusions: firstly, the scenarios, rightly, are different for different assets. Secondly, the scenarios depend on model assumptions beyond the seismic source models. In the application here, the model of the ground motion variability has a distinct effect on the scenarios for the residential-building stock. Given that the employed state-of-the-art model might overestimate the event-to-event variability (Bodenmann et al., 2023), the results for the building stock should be utilized carefully. Thirdly, because the earthquake scenario in the building case is representative of certain loss return periods only in combination with high ground motions, it should be investigated if and how the representative scenario should also provide ground motion fields together with the earthquake event. We note that these issues are also present for the classical hazard and loss disaggregation methods. Overall, for practical risk management tasks it is recommended to use the historic events jointly with the identified scenarios, in particular for the residential-building stock case.

The dependence of the results on the engineering system and the model assumptions implies that the representative scenarios should be regularly updated, depending on how much the analyzed system changes and the models improve. According to demographic projections of the Chilean National Statistics Institute (INE), based on the 2017 Chilean census, the population in the two communes will increase by around 13 % by 2035 with respect to the population in 2017 (INE, 2017). This demographical change will likely cause changes in the residential-building stock and the power demand. Depending on how these changes develop, the representative earthquake scenarios may change as well. Furthermore, the presented results do not consider the potential impact caused by further earthquake-triggered hazards, such as tsunamis and landslides. The tsunami impact during offshore earthquake scenarios with a large magnitude should be considered in a complete loss estimation and may affect the scenario selection associated with large return periods. Similarly, the scenario selection may change if one considers the landslides potential in the study area. Finally, cascading effects have not been considered in the power network model. Although the topology of the SEN and the redundancy of generators along the network reduce the probability of large blackouts, hub nodes far from the two communes may affect them through a sequence of failures triggered by earthquake scenarios impacting them. This may shift the representative scenarios to locations further from the study area.

7.2 On the definition of representative scenarios and the methodology for identifying them

We presented the evaluation of representative earthquake scenarios based on the loss occurrence and the loss exceedance approach; the latter coincides with the classical loss disaggregation method (Goda and Hong, 2009; Jayaram and Baker, 2009b). In the illustrative example of Sect. 4.2, we compare the results of the two approaches. For the case of hazard disaggregation, it has been proposed in the literature that the results of both approaches should be reported (Fox, 2023). However, we decided against reporting the scenarios of the exceedance approach for the communes of Valparaíso and Viña del Mar to avoid confusion. We find the loss occurrence approach to have a more intuitive interpretation. Scenarios identified with this approach correspond to a loss that is the t-year loss, which can be reported jointly with the scenarios. They are the most likely scenarios leading to this value (which on average is exceeded once in t years). By contrast, we find it difficult to communicate the meaning of the scenarios with the loss exceedance approach – and we believe it will be mostly misunderstood. Scenarios obtained with the loss occurrence approach can be described as “representative of a loss that is exceeded on average once in t years”. For the loss exceedance approach, one would need to describe scenarios as “representative of the losses that would occur when conditioning on a loss at least as large as the one that would be exceeded once in t years”, which seems too convoluted to communicate effectively. We also have difficulties in conceiving of a risk management activity for which such a definition would be more appropriate. However, we acknowledge that this discussion could benefit from additional comparisons of the two approaches in future studies.

To evaluate the representative scenarios, we adapted the methodology of Rosero-Velásquez and Straub (2022). The methodology leads to a lower computational cost in terms of loss evaluations compared to the classical loss disaggregation. By incorporating active learning, the methodology concentrates the conditional-loss evaluations around the scenarios that most likely produce the t-year loss value l_t. This concentration of samples around the solution and the smooth approximation of the conditional density with KDE make the methodology more suitable for selecting representative scenarios with the loss occurrence approach. For this approach, the classical loss disaggregation has to rely on the numerical derivative of the empirical CDF (Baker et al., 2021).

Although single representative scenarios are valuable for risk mitigation and communication purposes, they also have several limitations. For example, designing effective risk mitigation strategies, such as resource allocation before the event, using a single representative scenario would result in solutions tailored to the spatial distribution of damage of the specific selected scenario. Thus, better strategies could be defined by considering multiple scenarios, even for the same loss return period.

Possible extensions of the methodology include catalogs with multiple hazards (e.g., seismic scenarios with a tsunami), loss calculations considering indirect consequences, and high-dimensional scenarios (e.g., including the damage states of the individual components, either buildings or power network components). For the later, however, the dimensionality of the damage states has to be reduced (Rosero-Velásquez and Straub, 2019).

8 Conclusion

We present a methodology and algorithm to determine representative earthquake scenarios from a synthetic earthquake catalog. We applied the methodology to the communes of Valparaíso and Viña del Mar in Chile. Because the identified scenarios should be representative of extreme losses, they differ depending on the exposed assets. In this contribution, we consider the building stock and the electrical supply network. The application shows that the methodology can work and allows for the identification of scenarios more systematically than by selection or extrapolation from past events. However, the results for the residential-building stock also show that resulting scenarios cannot be considered independently of the resulting ground motions. Therefore, future work should investigate scenarios that also include the ground motions. Because the description of ground motion fields requires a large number of parameters, the existing methodology will need to be extended to be able to cope with such scenarios.

Appendix A: List of acronyms used in this paper

Acronym	Name
PGA	Peak ground acceleration
PDF	Probability density function
PSHA	Probabilistic seismic hazard analysis
GMM	Ground motion model
CDF	Cumulative distribution function
KDE	Kernel density estimation
GPR	Gaussian process regression
AL	Active learning
AEI	Augmented expected improvement
INE	Chilean National Statistics Institute
UNESCO	United Nations Educational, Scientific andCultural Organization
SEN	Chilean National Electric System
ONEMI	Chilean National Office of Emergency ofthe Ministry of the Interior and Public Security
USGS	United States Geological Survey
IS	Importance sampling
SARA	Project “South America Risk Assessment”
ENS	Energy not supplied
DCOPF	Direct current optimal power flow
USD	United States dollar
RP	Return period
ODEPLAN	Chilean National Planning Office
CSN	National Seismological Center of theUniversity of Chile

Appendix B: Variables used in this paper

Table B1 List of random variables used in this paper, distinguishing the random variable notation from a value and a sample value.

Table B2 List of indices used in this paper.

Table B3 List of numerical parameters used in this paper.

Table B4 List of variables and functions used in this paper. Variables with a hat (e.g., ${\widehat{l}}_{\mathrm{t}}$) represent approximations of their respective variables without a hat (e.g., l_t). Variables with the superscript exc refer to the exceedance approach, and the corresponding variables without this superscript correspond to the occurrence approach.

Table B5 List of additional variables used in this paper for the examples and applications in Sects. 4 and 5.