# Ground motion simulations in Quito (Ecuador) due to major earthquakes from the subduction zone

Françoise Courboulex, David Alejandro Castro Cruz, Aurore Laurendeau, Luis Fabian Bonilla, Alexandra Alvarado, Etienne Bertrand

► To cite this version:

Françoise Courboulex, David Alejandro Castro Cruz, Aurore Laurendeau, Luis Fabian Bonilla, Alexandra Alvarado, et al. Ground motion simulations in Quito (Ecuador) due to major earthquakes from the subduction zone. Geophysical Journal International, 2022, 229 (3), pp.2192-2208. 10.1093/gji/ggac044. hal-03617258

HAL Id: hal-03617258

https://hal.science/hal-03617258

Submitted on 25 Jun 2022

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
<figure><img src="image_2.png" /></figure>

# Ground motion simulations in Quito (Ecuador) due to major earthquakes from the subduction zone

Francoise Courboulex<sup>1</sup>,<sup>5</sup> David Alejandro Castro-Cruz<sup>1,2,†</sup> Aurore Laurendeau<sup>1,5</sup> Luis Fabian Bonilla<sup>4</sup>,<sup>3</sup> Alexandra Alvarado<sup>3</sup> and Etienne Bertrand<sup>2,4</sup>

Accepted 2022 February 2. Received 2021 November 17; in original form 2021 April 28

## SUMMARY

In 1906, an earthquake with a magnitude estimated between $M_w$ 8.4 and 8.8 occurred in the subduction zone along the coast of Ecuador and Colombia. This earthquake caused extensive damage on the coast but had a rather small impact on the capital city of Quito, situated 180 km away. At that time, the city of Quito extended over a small area with a few thousand inhabitants, while today it stretches over 40 km and has a population of over 3 million, with most of the city built without paraseismic regulations. The aim of this study is to obtain new insights on the impact that large earthquakes from the subduction zone would have on the city today. This question is crucial since we know that the city of Quito is prone to site effects and that the southern part of the city amplifies seismic waves at low frequencies, around 0.3–0.4 Hz. In April 2016, an $M_w$ 7.8 earthquake occurred on the subduction interface in the Pedernales area. This event was the first large earthquake in the city of Quito to be well recorded by 13 stations of the permanent accelerometric network (RENAC). In this study, we take advantage of this data set (main shock and large aftershock recordings) to (1) test an empirical Green’s function blind simulation approach where the input stress drop is taken from a global catalogue of source time functions, (2) compare the synthetic accelerograms and ground motion values we obtain for an $M_w$ 7.8 earthquake with the actual recordings of the Pedernales earthquake and then (3) simulate larger earthquakes of $M_w$ 8.2 and 8.5 from the subduction zone. For $M_w$ 7.8 simulations, our approach allows a good reproduction of the ground motions in the whole frequency bands and properly takes into account site effects. For $M_w$ 8.2 and 8.5 simulations, we obtain for the stations in the southern part of the basin, larger values at low frequencies than the predicted motion given by ground motion models. These values, although high, should be supported by new or recent buildings if they are constructed respecting the building code that applies in Quito. Therefore, for this type of strong but distant earthquake, the seismic standards appear to be well suited and it is imperative to ensure that they are well considered in the design of the new buildings to be constructed, especially in the southern part of the expanding city.

**Key words:** South America; Earthquake ground motions; Earthquake hazard; Site effects.

## 1 INTRODUCTION

The country of Ecuador is prone to a high seismic hazard that mainly results from the subduction of the Nazca oceanic plate under the South American continental plate at a rate of ~6 cm yr⁻<fnref n="1" />

(Kendrick et al. 2003). Since the beginning of the 20th century, five earthquakes greater than $M_w$ 7.5 have occurred along the Ecuador–Colombia subduction zone (Font et al. 2013; Nocquet et al. 2014, 2017; Yepes et al. 2016). Their epicentres are represented in Fig. 1 by red circles. The largest one was the 1906 Esmeraldas mega-earthquake, whose magnitude was evaluated at $M_w$ 8.8 by Kanamori & McNally (1982), and $M_w$ 8.4 by a recent analysis of tsunami data

<footnote n="1">¹Université Côte d'Azur, CNRS, IRD, Observatoire de la Côte d'Azur, Géoazur, 06560 Valbonne, France. E-mail: courboulex@geoazur.unice.fr</footnote>

<footnote n="2">²CEREMA, Agence de Sophia Antipolis, 06903 Valbonne, France</footnote>

<footnote n="3">³Instituto Geofísico, Escuela Politécnica, Ladrón de Guevara E11-253, Apartado 2759, Quito, Ecuador</footnote>

<footnote n="4">⁴GERS Department, Université Gustave Eiffel, 77447 Marne-la-Vallée, France</footnote>

<footnote n="5">⁵Institut de Radioprotection et de Sûreté Nucléaire (IRSN), PSE-ENV, SCAN, BERSSIN, 92262 Fontenay-aux-Roses, France</footnote>

<footnote>†Now at: King Abdullah University of Science and Technology, Saudi Arabia.</footnote>
<figure><img src="image_3.png" /><figcaption>**Figure 1.** Ecuador main historical seismicity (Beauval *et al.* 2013) and tectonic settings. Red circles: epicentres of superficial earthquakes (depth <50 km.) with $M_w \ge 7.0$ for subduction earthquakes and $M_w \ge 6$ for crustal earthquakes. Q indicates the location of the capital city of Ecuador, Quito. The light-red ellipsoids on the coast represent the approximate rupture area of the earthquakes simulated in this paper. The smallest correspond to the $M_w 7.8$ events and the largest to the $M_w 8.5$. They are there only to give an idea of the area that can be broken by the earthquakes, but this size strongly depends on the rupture parameters. The two stars represent the position of the earthquakes used as EGFs for the simulation of $M_w < 8$ events. Figure modified from Alvarado *et al.* (2014).</figcaption></figure>

(Yoshimoto *et al.* 2017). The last large earthquake occurred on 2016 April 16 (Pedernales earthquake, $M_w 7.8$). It broke the interplate interface at shallow depth along the subduction zone and caused extensive damage on the coast (Ye *et al.* 2016; Nocquet *et al.* 2017). It was the first large earthquake in Ecuador to be well recorded by seismic stations throughout the country, including in the capital, Quito.

Quito is a crowded city (3 million inhabitants) situated in a narrow valley (~5 to 8 km wide and 40 km long) perched in the Andes Cordillera (2400–3000 m altitude) and 180 km from the coast (letter Q in Fig. 1). The known historical earthquakes that have damaged Quito were located on the faults of the Cordillera, such as the Guaylabamba 1587, Riobamba 1797, Quito 1859 and Ibarra 1868 earthquakes (Del Pino & Yepes 1990; Beauval *et al.* 2013; Alvarado *et al.* 2014). The city is built on the hanging wall of an active reverse fault that generates moderate size earthquakes (Alvarado *et al.* 2014; Vaca *et al.* 2019) and is partially creeping (Marinière *et al.* 2020). The seismic hazard related to the activity of this fault should not be neglected, but it is not the topic of this study.

Our aim is to quantify the ground motions that could be generated in Quito by large earthquakes ($M_w \ge 8$), occurring in the subduction zone. When the Esmeraldas earthquake occurred in 1906, Quito was a city with a small extension that cannot be compared with the present day agglomeration of 3 million inhabitants (see the evolution of the city at https://comet.nerc.ac.uk/tomorrows-cities-quito/). Quito is situated rather far away from the subduction interface (180 km), but we know that this distance is not large enough to ensure that the city is safe against a strong earthquake from the subduction zone. Indeed, Mexico City suffered heavy damage in 1985 after an $M 8.1$ event that occurred in the subduction zone 400 km away because seismic waves were strongly amplified by the sedimentary basin underlying the city (Anderson *et al.* 1986). Are the lithological site conditions also capable of amplifying seismic waves in the Quito basin?

Quito is built on a piggy-back basin filled with volcanic and fluvial origin deposits, the structure of which is still poorly known. The actual knowledge of the basin is limited to a few tens of metres, thanks to stratigraphic outcrop observations (Alvarado *et al.* 2014)
and some geotechnical and geophysical investigations. Laurendeau et al. (2017) analysed the recordings of ~200 earthquakes at 18 stations of the permanent accelerometric network RENAC (Red Nacional de Acelerógrafos, operated by Instituto Geofísico de la Escuela Politécnica Nacional) and found that the southern part of the basin strongly amplifies seismic waves at frequencies around 0.35 Hz (by a factor of 3–5). At higher frequencies ($f > 1$ Hz), they also obtained significant amplifications but with a large spatial variability, coherent with the results obtained by ambient vibration analysis (Chatelain et al. 1999; Guéguen et al. 2000; Reyes Once 2020). Site amplification, although still poorly understood, is now well established (Pacheco et al. 2022). It is therefore important to evaluate its impact when a future large earthquake occurs.

The easiest way to predict ground motion due to future earthquakes is to select appropriate ground motion models (GMMs) with ad hoc parameters for source and site effects. This approach is widely used by seismologists and earthquake engineers to assess the seismic hazard in earthquake prone regions (e.g. Beauval et al. 2018). However, GMMs cannot properly take into account particular site effects, and they are not able to produce synthetic seismograms, which is of the utmost importance in cities like Quito, where few or no large earthquakes have yet been recorded. Ground motion simulations can then be undertaken to produce synthetic seismograms, which can be used in earthquake scenarios, PSHA (probabilistic seismic hazard assessment), or the modelling of the dynamic behaviour of buildings. Several methods exist, each with its advantages and disadvantages (e.g. Douglas & Aochi 2008 and references therein). The recent growth in computational power now makes possible the numerical simulation of broad-band seismograms that take into account the latest knowledge of the source process and a refined description of the propagation media at large and small scales (e.g. Cruz-Atienza et al. 2016; Gallovič 2017; Rodgers et al. 2018; Ulrich et al. 2019).

However, there are few places in the world where the ground shallow structure has been sufficiently well characterized to perform such numerical studies. Another issue is that numerical simulations often struggle to constrain the variability of their results when numerous input parameters are needed (i.e. slip distribution, rupture velocity, stress drop, rise time, anelastic attenuation, site effects, etc.).

Pulido et al. (2020) proposed for the first time a broad-band simulation of the 1906 mega-earthquake. They converted the results of their simulations into macroseismic intensities and found a rather good agreement with the values reported from all over the country for the 1906 earthquake. They considered a single point for Quito city and did not aim to refine the effects due to site effects in different parts of the city.

In order to take into account site effects in Quito and to reduce the variability due to numerous input parameters, we use a well-constrained empirical Green's function summation method. In this method, the recordings of relatively small earthquakes are used to account for path and site effects in a realistic way. The main input parameter to constrain, the stress drop of the target earthquake, is chosen from the source time function (STF) duration of the global SCARDEC database (Vallée et al. 2011; Vallée 2013; Vallée & Douet 2016). The method is presented in Section 2 and the data used in Section 3.

In a first step, we simulate the ground motions in Quito generated by several $M_w$ 7.8 earthquakes from the subduction zone (with variable STF duration). We locate the target earthquakes at the place where the 2016 Pedernales event occurred, in order to be able to compare the simulation results with actual recordings of the

Pedernales earthquake. We select four aftershock recordings ($M_w$ 5.8–6.1) as empirical Green's functions (EGFs). The simulation results and comparison with the Pedernales earthquake recordings are analysed in Section 4.

In a last step (Section 5), we simulate the ground motion generated by larger earthquakes ($M_w$ 8.2 and 8.5) in order to highlight the main ground motion values expected in the city of Quito when such earthquakes will occur.

## 2 METHOD

### 2.1 Empirical Green's function method

Since the geometry and velocity of the sediments are poorly known in the Quito basin, as well as the regional velocity model between the subduction zone and Quito, we chose to use the recordings of small earthquakes as EGFs. This approach takes directly into account both path and site effects at different stations, as they are included in the EGFs (Hartzell 1978). We assume that the recordings of the small earthquakes represent the Green's functions at each point of the fault plane activated during the rupture of the large simulated event (called the target earthquake).

We use a stochastic EGF method that requires few input parameters and a minimal description of the source. The method was proposed by Kohrs-Sansorny et al. (2005) following the work of Wennerberg (1990), Ordaz et al. (1995) and Joyner & Boore (1986). Since then, it has been used and improved by several authors (e.g. Courboulex et al. 2010; Salichon et al. 2010; Honoré et al. 2011; Vieux-Champagne et al. 2017; Wang et al. 2017, 2019; Lorenzo et al. 2018; Niño et al. 2018) to simulate the ground motions due to future earthquakes in different environments.

The method consists in generating a large number of synthetic STFs (called equivalent STFs, eSTFs) and to convolve them with the EGFs at each station. Each eSTF represents the specific rupture of a given earthquake up to the corner frequency of the EGF. Only four parameters are needed to build the eSTFs: the moment and corner frequency of the small earthquake used as an EGF ($m_o$ and $f_c$) and of the target earthquake ($M_o$ and $F_c$). The eSTFs are built using a random process based on probability densities as proposed by Ordaz et al. (1995), which ensures the generation of the $\omega^{-2}$ model spectra at all frequencies (Brune 1970). The source frequency content must, on average, explain the theoretical ratio $R(f)$ between the spectra of the large and small events (see Fig. A1 of the Annex section). If the two earthquakes follow the $\omega^{-2}$ model then

$$R(f) = \frac{M_o}{m_o} \cdot \frac{1 + \left(\frac{f}{f_c}\right)^2}{1 + \left(\frac{f}{f_c}\right)^2} \quad (1)$$

The generation of the eSTFs is performed in two stages as proposed by Khors-Sansorny et al. (2005) to improve its variability: in the first stage, a number $\eta_c$ of delays $t_c$ are randomly generated with a probability density $\rho(t)$ over the whole source duration $T$. In the second stage, a number $\eta_d$ of delays $t_d$ are again generated with a second probability density over a window duration $T_d \le T$ centred on each delay generated in the first stage. Finally, $\eta = \eta_c \cdot \eta_d$ small events are summed together and scaled by a factor $\kappa$. The choice of the parameters $\eta$ and $\kappa$, the probability density and all the details in the generation of the eSTF are given in the Annex section.
## 2.2 Stress drop variability constrained by a global STF database

The stress drop ratio between the large and the small earthquake can be expressed as $C = M_0 F_c^3 / m_0 f_c^3$. C strongly influences the level of acceleration obtained at the surface as shown by Honoré et al. (2011) and analysed by Wang et al. (2019). It is then important to correctly estimate its value and variability. For a given EGF, C directly depends on the value of $F_c$, which is also directly linked, through a rupture model, to the duration T of the STF.

We propose to determine the mean value and the variability of T by using the SCARDEC global database of STFs (Vallée et al. 2011, 2013; Vallée & Douet 2016). SCARDEC is a consistent database of source parameters (seismic moment, focal mechanism and STF) of earthquakes of magnitude larger than 5.8 since 1992. Courboulex et al. (2016) analysed the SCARDEC database and showed that the duration of the STF is a log-normally distributed quantity for a given magnitude. The log of its value can then be represented by a Gaussian curve with a mean value and a standard deviation that represents its variability. For a given value of $M_w$, Courboulex et al. (2016) and Chounet & Vallée (2018) found that the STF duration is larger for earthquakes in the subduction zone than outside (Fig. 2). A simple regression makes it possible to obtain the following mean duration T for earthquakes in and outside the subduction zone (Courboulex et al. 2016):

$$ \log 10(T) = 0.2867 \log 10(M_0) - 4.3205 (\text{SUB}), \quad (2) $$

$$ \log 10(T) = 0.3101 \log 10(M_0) - 4.9032 (\text{NOT-SUB}). \quad (3) $$

The standard deviation (hereafter called sigma) values are 0.32 for subduction events and 0.34 for others in natural log (respectively 0.138, 0.147 in log10), without any dependence on the magnitude. These values are smaller than the variability usually reported for stress drop values in global databases (Cotton et al. 2013; Courboulex et al. 2016).

We use these results (eq. 2 for subduction zone and corresponding sigma value) in our procedure (Fig. 3) to build a large number of eSTFs that account for the natural variability of the rupture processes. The details of the procedure are given in the Annex section and the corresponding python code is available (see data availability). The eSTFs are then convolved with the EGFs to obtain a large number of synthetic seismograms at each station and each component. From those, we compute PGA, Fourier and response spectra.

# 3 DATA

## 3.1 The 2016 M<sub>w</sub> 7.8 Pedernales earthquake

The largest recent earthquake in Ecuador (Pedernales earthquake, 2016 April 16, M<sub>w</sub> 7.8) caused extensive damage to cities on the coast and the death of 700 persons. This event occurred on the subduction interface at shallow depth on a part of the fault that already broke during the 1942 earthquake.

The rupture area has been well identified thanks to seismologic, geodetic and tsunami data (Nocquet et al. 2016; Ye et al. 2016). It began at the north and then propagated almost unilaterally towards the south with an area of higher slip located 70 km south of the hypocentre (Nocquet et al. 2016). Many studies have been undertaken on this earthquake in order to localize the aftershocks (Agurto-Detzel et al. 2019; Meltzer et al. 2019), specify the structure (León-Ríos et al. 2019) and understand the relationships between coseismic slip, slow

slip, afterslip and microseismicity (Rolandone et al. 2018; Tsang et al. 2019; Chalumeau et al. 2021). Because of the significant distance (~180 km) from Quito to the source, it is the low frequencies of the signal that were felt, producing oscillations of buildings without serious apparent damage.

It is the first large event to be recorded on all the RENAC seismic networks of the country (Beauval et al. 2017). Laurendeau et al. (2017) analysed the recordings of the main shock in Quito and found a clear amplification of the seismic waves around a frequency of 0.3–0.4 Hz, similar to the values they already obtained for smaller earthquakes. This amplification is obvious on accelerograms (see fig. 5 of Laurendeau et al. 2017) and can be a subject of preoccupation in the case of a larger earthquake.

We use the recordings of the Pedernales earthquake and its main aftershocks on the permanent accelerometric stations of Quito to test the blind simulation methodology described above. Since the city is located 180 km from the coast and the rupture is almost parallel to the main elongation of the city (Fig. 4), we have a simple configuration that avoids the influence of directivity effects of the rupture process on the ground motions. The large distance also makes it possible to neglect near source effects and thus justifies the use of a simple method.

## 3.2 Choice of the aftershocks used as EGFs

The Pedernales earthquake was followed by a large number of aftershocks (Agurto-Detzel et al. 2019; Meltzer et al. 2019). Because of the large epicentral distance, the low sensitivity of the accelerometric stations in Quito, and the high level of noise produced by anthropic activity, only the largest aftershocks ($M_w > 5.7$) were well recorded with a high signal-to-noise ratio. The frequency band we select is 0.1–10 Hz for all the study. At larger frequencies, the city noise is much too important and at lower frequencies there is almost no signal recorded by the stations. We have kept four very well-recorded aftershocks with a signal-to-noise ratio larger than 5 in all the considered frequency bands: two were located to the north (EGF-N1 and EGF-N2) of the Pedernales rupture and two were located to the south (EGF-S1 and EGF-S2) (see Fig. 4 and Table 1). We use the recordings of each aftershock successively as EGFs, in order to estimate the variability induced by this choice. We apply an attenuation correction (geometric attenuation only) to move the aftershock locations to the place where the Pedernales earthquake released the larger part of its energy (see the dotted ellipse on the left-hand part of Fig. 4). For each earthquake, the corner frequency is determined carefully by fitting with a theoretical $\omega^{-2}$ model (Brune 1970) in displacement and acceleration (the characteristics of the four earthquakes taken as aftershocks are in Table 1).

# 4 SIMULATION OF M<sub>w</sub> 7.8 EARTHQUAKES FROM THE SUBDUCTION ZONE

## 4.1 Target event characteristics

As for a blind simulation, we only impose the magnitude and the location of the centroid of the target event. We choose a magnitude of M<sub>w</sub> 7.8 and the location at the place where the main Pedernales slip occurred on the fault (see ‘maximum slip zone’ in Fig. 4) to be able to compare the simulations with the actual Pedernales
<figure><img src="image_4.png" /><figcaption>**Figure 2.** Left: source time function duration versus magnitude ($M_w$) for earthquakes that occurred on the subduction interface (grey squares, SUB) and others (black diamond, NOT-SUB) from Courboulex *et al.* (2016). The red square represents the STF duration of the 2016, $M_w$ 7.8 Pedernales earthquake. Right: STF of the Pedernales earthquake from SCARDEC database (Vallée & Douet 2016).</figcaption></figure>

<figure><img src="image_5.png" /><figcaption>**Figure 3.** Schematic view of the procedure. A large number (here 500) of eSTF are generated using three input parameters (magnitude of the target event $M_w$, magnitude $m_w$ and corner frequency $f_c$ of the smaller event used as the EGF) and the STF duration variability given by the SCARDEC data base. The eSTF generation follows the two-step procedure of Kohrs-Sansorny *et al.* (2005). The eSTFs are convolved with the EGFs at each station to obtain the synthetic seismograms.</figcaption></figure>

recordings. Nothing is imposed concerning the stress drop, rupture velocity, rupture length or maximum slip.

In the process, we choose eq. (2), for subduction zone earthquakes, to generate a large number of equivalent eSTF. Fig. 5 (top) shows the log-normal distribution of the duration of 500 eSTF. In a few cases, the source duration is very short (22 s), and in a few cases it is very long (80 s); most of the eSTFs have a duration between 30 and 60 s, following what is obtained on the STF SCARDEC database (Courboulex *et al.* 2016) for an $M_w$ 7.8 earthquake. We also show in Fig. 5 the corresponding stress drop values using the relation $\Delta\sigma = 7/16 M_0 F_c^3 / (k Vs)^3$, with $k = 0.23$ and $Vs = 3900$ m s$^{-1}$. Vs is the shear wave velocity near the source and k a parameter introduced by Madariaga (1976) that depends on the rupture model and the type of wave. $F_c$ is taken equal to $0.8/T$ ($T$ being the STF duration). It is important to remember that absolute values of the stress drop strongly depend on the parameters used in the equation and should not be further interpreted (see fig. 3 of Courboulex *et al.* 2016).

## 4.2 Results and comparison with the Pedernales earthquake recordings

Each eSTF is convolved with the recordings of the small earthquake selected as an EGF, at each station and each component, in order to produce a large number of synthetic accelerograms that can be treated statistically. Fig. 6 presents an example of a few synthetic accelerograms obtained with different eSTF durations at four stations in the Quito basin using EGF-S1. Station 24MA is considered by Laurendeau *et al.* (2017) as the reference station with no site effects (its H/V ratio is almost flat), ALLO is situated outside the Quito basin at a high elevation point, CIRC is in the north part of the basin and LILI in the south.

For a given station, the large variability obtained in the waveforms (Fig. 6) is due to the variability in eSTF duration. The source process of the Pedernales earthquake had a total duration of 43 s (see Fig. 2, right side). In Fig. 6, we see that it is well consistent with the synthetic on line 2 or 3, for which eSTF has a similar duration.
<figure><img src="image_6.png" /><figcaption>**Figure 4.** Left: map of Ecuador and the Pacific Ocean (the highest altitude zones are in dark grey and the ocean is in white). The projection of the slip zone of the Pedernales $M_w$ 7.8 earthquakes from Nocquet *et al.* (2017) is coloured. Black stars indicate the location of the four aftershocks used as EGFs. Right: zoom on the accelerometric stations (RENAC network) of the Quito basin used in this study.</figcaption></figure>

<table><caption>**Table 1.** Characteristics of the aftershocks selected as EGFs. Location, initial time and $M_w$ are from the USGS catalogue. EGF-N1 and EGF-N2 are located to the north of the Pedernales rupture zone. EGF-S1 and EGF-S2 are located to the south (Fig. 4).</caption><thead><tr><th rowspan="2">M<sub>w</sub></th><th rowspan="2">EGF Name</th><th colspan="3">Origin time (UTC)</th><th rowspan="2">Location Longitude (°)</th><th rowspan="2">Depth (km)</th><th rowspan="2">F<sub>c</sub> (Hz)</th></tr><tr><th>Date</th><th>Hour</th><th>Latitude (°)</th></tr></thead><tbody><tr><td><b>5.7</b></td><td>EGF-N1</td><td>19/04/2016</td><td>22:22</td><td>0.578</td><td>-80.025</td><td>15.4</td><td>1.15E-01</td></tr><tr><td><b>6.0</b></td><td>EGF-S1</td><td>22/04/2016</td><td>3:03</td><td>-0.292</td><td>-80.504</td><td>10.0</td><td>1.06E-01</td></tr><tr><td><b>5.9</b></td><td>EGF-S2</td><td>22/04/2016</td><td>3:25</td><td>-0.281</td><td>-80.504</td><td>10.3</td><td>9.00E-02</td></tr><tr><td><b>5.7</b></td><td>EGF-N2</td><td>23/04/2016</td><td>1:24</td><td>0.613</td><td>-80.252</td><td>10.0</td><td>1.10E-01</td></tr></tbody></table>

The process is applied using the four EGFs successively at each component of each station. There are numerous synthetic traces (500 eSTF, 13 stations, 3 components, 4 EGFs, i.e. a total of 78 000 traces). We analyse the results in time, in the Fourier domain and on the elastic response spectra. This last representation has the advantage of giving an integrated vision of the results in different frequency bands and also to give an idea of the impact that the signals would have on buildings.

From each synthetic accelerogram, we extract the peak ground acceleration (PGA) and represent it in function of the eSTF duration (Fig. 5, for the case of station LILI). Note that the accelerograms are filtered between 0.1 and 10 Hz, so this PGA value does not integrate higher frequencies. The distribution of PGA is shown on the right-hand part of Fig. 5. PGA value decreases when the eSTF duration becomes larger, following the equation:

$$ \log(\text{PGA}) = a \log(T) + b, \quad \text{with } a = -2.20 \pm 0.03 \text{ and } b = 3.2 \pm 0.1 \qquad (4) $$

The value of the slope, *a*, is almost similar to the proportionality relation: $\log(\text{PGA}) \sim 2.4 \log(f_c)$, obtained by Causse & Song (2015). As expected, when the total rupture duration is smaller (or

$f_c$ higher), the energy is concentrated in a shorter time and then the level of acceleration is higher.

From Fig. 5, we see that the PGA variability is mainly due to the eSTF duration variability. The choice of the EGF also has an impact on the results but it is rather low and is not related to the location of the EGFs to the south or north of the rupture. This independence is certainly due to the special situation of Quito (see Fig. 1), parallel to the subduction zone and rather far away from it. The PGA variability given by the standard deviation, sigma, of the log-normal distribution of the PGA is 0.62 (natural log) for this station. This value is comparable to the value given by GMMs (from 0.5 to 0.8; see Cotton *et al.* 2013).

The same simulations are represented in the Fourier spectral-domain (Fig. 7). For clarity, only the mean (continuous lines) and standard deviations (dashed lines) are represented. Interestingly, the specific shape of the Fourier spectra is well reproduced by the simulations at each station (Fig. 7). We observe some differences at frequencies around 0.13 Hz, especially for the simulation with EGF-N2 that cannot be explained simply. This difference is certainly representative of the variability that exists in the natural process.

The same simulations are represented (Fig. 8) in terms of the elastic response spectrum (considering 5 per cent damping). The
<figure><img src="image_7.png" /><figcaption>**Figure 5.** Bottom: example of one stochastic distribution of eSTF duration *T* (s) that follows on average a log-normal distribution (this corresponds to the 'random generation of *T* values' step in the procedure described in the diagram of Fig. 3). Top left: PGA values of the synthetic signals at station LILI (NS component) in function of the eSTF duration, for the four different EGFs. The corresponding stress drop values Δ*σ* are indicated with *k* = 0.23 (Madariaga 1976) and Vs = 3900 m s⁻¹. The distribution of eSTF duration (bottom) and PGA (right) are presented with the mean and sigma values depicted by red lines.</figcaption></figure>

spectral shapes at the 24MA and CIRC stations are rather similar but it is clearly different for ALLO, exhibiting a high narrow peak around a period of 0.5 s, and LILI with a broad-band frequency spectra and significant values at a period around 2 to 3 s, much higher than at the other three stations. These features, due to local site effect, are also present in the Pedernales recordings (bold orange line). For all the stations and all the EGFs, the Pedernales recordings fall in the range predicted by the mean ± one standard deviation predicted by the simulations.

For comparison, we also plot in Fig. 8 (grey shaded zone) the prediction of the Abrahamson *et al.* (2016) GMMs (hereafter called AB2016) developed for interface subduction earthquakes in the specific case of the backarc region. Indeed, Beauval *et al.* (2017) found it particularly adapted to predict the 2016 Pedernales earthquake recordings and two large aftershocks. The Vs30 parameter is fixed equal to an average value of 400 m s⁻¹ for each station. We see that, for stations CIRC and 24MA, the AB2016 prediction fits rather well the actual Pedernales data and our simulations. For stations ALLO and LILI, which are prone to site effects, we see that our simulations better fit the Pedernales observations than AB2016. This is not surprising as it well known that GMMs may not adequately capture site amplification effects.

The results obtained at each station are now represented on Fig. 9 for three selected periods: *T* = 0.1 s (equivalent to the PGA at 10 Hz), *T* = 0.6 s, the period that gives the highest value of the response spectrum on most stations, and the period *T* = 2.8 s, which corresponds to the second peak visible on LILI (and the other stations in the southern part of the basin). For comparison, we also plot the lines that correspond to the mean and mean ± sigma values of AB2016. The distances used for the equation and indicated in Fig. 9 are the closest distance to the rupture. First, in Fig. 9 we observe that, for the three natural periods, the real values of the

Pedernales earthquake (orange star) always fall within the error bars predicted by our simulations, and this, whatever the chosen EGF. This indicates that the choice of the EGF, in this configuration, does not have a strong impact on the ground motions levels.

For high frequencies (PGA, Fig. 9, top), although the values are quite variable from one station to another, our simulations (mean ± sigma) always include the values recorded during the Pedernales earthquake. There is one exception for station SADP, which is far north outside the basin. The AB2016 prediction also fits rather well the Pedernales PGA values, but it slightly overestimates the values for stations PRAM, AEPN, CMEA CIRC and 24MA located in the middle or north of the basin.

At a period *T* = 0.6 s, our simulations almost perfectly reproduce all the Pedernales values. AB2016 predictions are very good for the stations PRAM, AEPN, CMEA CIRC and 24MA and still good for the other stations where the Pedernales values remain within the confidence interval of the predictions.

For a period $T = 2.8$ s, the recorded values of the Pedernales earthquake at four stations in the south QUIB, ZALD LILI and VILF are significantly higher than the value predicted by AB2016 (and it is not the Vs30 value that can change anything at this frequency). Our simulations reproduce this low-frequency remarkably well.

In general, we always obtain an excellent agreement, in a large frequency band, between the results of our simulations and the actual data from the Pedernales earthquake whatever the chosen EGF. These simulations better fit the shape of the Pedernales earthquake than AB2016, especially for the stations situated in the southern part of the Quito basin. Moreover, the variability imposed by our method is not higher than that of GMMs (Fig. 9), which makes it useful in a blind simulation case. The next step is then to simulate stronger earthquakes for which we have no record.
<figure><img src="image_8.png" /><figcaption>**Figure 6.** Example of synthetic accelerograms ($M_w$ = 7.8) obtained on the NS component of four stations in the Quito basin (24MA, ALLO, CIRC and LILI). The value to the right of each trace indicates the duration of the eSTF used for this simulation. The real Pedernales accelerograms are represented at the bottom for comparison (orange line). It fits well with the synthetics of lines 3 and 4, whose eSTF durations are similar to the Pedernales source duration (43 s). The EGF-S1 derived from the earthquake of 2016 April 22 was used.</figcaption></figure>

## 5 SIMULATION OF $M_w > 8$ EARTHQUAKES ALONG THE SUBDUCTION ZONE

As already mentioned, the $M_w$ 7.8 Pedernales earthquake caused severe destruction in the country but only slightly damaged the city of Quito. An important question is whether a stronger earthquake from the subduction zone would threaten the capital or not. In this assessment, the amplification effect identified in the southern Quito City basin for low frequencies should not be neglected. It is therefore essential to continue to use empirical Green's functions to properly account for site effects.

In 1906, one mega-earthquake occurred along the subduction zone. It broke a fault which, according to Kelleher's work (1972) from macro-seismic intensities, has a total length under 500 km with the maximum slip to the north, at the level of Colombia, quite far from Quito. Its magnitude $M_w$ was calculated by Kanamori & McNally (1982) to be equal to 8.8. Ye et al. (2016) adopt a magnitude $M_w$ ~8.6 compatible with the $M_s$ = 8.6 value determined by Gutenberg & Richter (1954) with ancient seismograms and tsunami magnitude. Recently, Yamanaka et al. (2017) re-analysed the tsunami

signals of this earthquake and found $8.2 < M_w < 8.6$ with a dominant slip also in the northern part near the coast of Colombia.

Although Sedaghati et al. (2020) used the same type of point-source method to simulate the $M_w$ 9 Tohoku earthquake, we think that our approach is not valid for such large events. Indeed, earthquakes with a magnitude greater than 8.5 being rare in the SCARDEC database, it is rather difficult to constrain the mean duration of their STF. Moreover, our EGF method would be difficult to apply to a very large fault, since we do not have any recordings of earthquakes so far north that could be used as EGFs.

We choose to present the results of $M_w$ 8.2 (lowest magnitude estimated for the 1906 earthquake) and $M_w$ 8.5 (upper limit for the method and the source-city situation) earthquakes that would occur on the subduction interface (see dashed zones on Fig. 1). Because the impact of the choice of the EGF is rather low, we only present the results obtained with EGF-N1 and EGF-S1. In the first scenario (south-scenario), the main slip zone is located in the southern part of the Pedernales rupture area (southern star of Fig. 1), and then EGF-S1 is selected. In the second scenario (north-scenario) the main slip zone is located in the northern part of the Pedernales rupture (centre star on Fig. 1), and then EGF-N2 is selected. Of course the
<figure><img src="image_9.png" /><figcaption>**Figure 7.** Simulation results ($M_w$ 7.8) at four stations and four different EGFs represented as Fourier acceleration spectra (FAS) in accelerations (NS component). Mean values of the simulations are in a continuous line and mean value ± sigma are represented by dashed lines. The real recordings of the Pedernales earthquake are represented by an orange line. Frequencies lower than 0.1 Hz and higher than 10 Hz have been filtered and thus cannot be interpreted.</figcaption></figure>

waveform of an event further north (in the area where the $M_w$ 8.1 earthquake initiated in 1979, see Fig. 1) would have been interesting also, but there is still no event available from this zone that is well recorded on the Quito network.

The results are shown for all the stations for PGA and PSA values at periods of 0.1, 0.6 and 2.8 s (Fig. 10) and are also presented in the time domain (Fig. 11) and response spectra (Fig. 12) for four stations for $M_w$ 8.5. We observe first in Fig. 10 that the values obtained with the two scenarios are quite similar and that the main variability is still related to the STF duration. However, for low pseudo-periods (PGA and PSA at 0.1 s), the north-scenario generates slightly larger values on the stations in the northern part of Quito, like SADP, FENY, PRAM and CRPG, and the south-scenario generates larger values for the stations in the southern part. This tendency tends to disappear for higher pseudo periods. For PSA = 2.8, the north-scenario provides higher values for all the stations. This means that for this period, the effect of distance is minor and that the frequency content of the EGF is certainly predominant. For comparison, again, we plot on the same figure the AB2016 prediction.

For high frequencies (PGA and PSA at 0.1 s), our mean simulations provide values slightly lower than AB2016 for the $M_w$ 8.2 case and similar values for the $M_w$ 8.5 case. The mean PGA

values expected are between 0.5 and 1 m s⁻² for $M_w$ 8.2 and between 0.7 and 2 m s⁻² for $M_w$ 8.5. The largest values are obtained at station ALLO, which is prone to high-frequency site effect, and station SADP situated outside the Quito basin. These stations have resonance frequencies respectively at 2 and 1 Hz (see Laurendeau *et al.* 2017, supplement), which explains these high PGA values.

For a period $T = 0.6$ s, the simulations are generally higher than predicted by AB2016. This is particularly clear for $M_w$ 8.5, where the mean values of the simulations are close to the mean value + one sigma of AB2016. Again, the values are stronger at stations SADP and ALLO.

For a period of $T = 2.8$ s, also, the mean values predicted by the simulations are all higher than AB2016 for $M_w$ 8.2 and even more for $M_w$ 8.5. This tendency is particularly clear for the stations in the southern part of the basin (QUIB, ZALD, LILI and VILF) for which the simulations are higher than the mean + one sigma value of AB2016. The variability in the synthetic signals at each station due to different eSTF can be illustrated by the example shown in Fig. 11 on four stations for $M_w$ 8.5. The associated response spectra (Fig. 12) for the north and south scenarios are compared with AB2016 (grey shaded area). It is easy to see in this figure that for PGA, our simulations give the same predictions as AB2016, except
<figure><img src="image_10.png" /><figcaption>**Figure 8.** Same results as Fig. 7 represented as elastic response spectra versus period (damping 5 per cent). BC-hydro Abrahamson *et al.* (2016) GMMs' predictions (AB2016) are also presented (for Vs30 = 400 m s⁻¹) by grey zones (mean ± sigma).</figcaption></figure>

for station ALLO. This is the contrary for periods higher than 1 s: our simulations give larger values than AB216, but similar results for ALLO. The particular shape of the response spectra obtained at station LILI (similar to stations VILF, QUIB and ZALD) clearly show the difference with the AB2016 prediction.

# 6 DISCUSSION

This paper presents the first attempt to simulate the ground motion that would be generated in the city of Quito by large earthquakes ($M_w > 8$) from the subduction interface. This simulation is motivated by the mega-earthquake that took place in 1906 along the Colombian–Ecuadorian coast, but does not attempt to reproduce its effects, because at that time the city of Quito was small and had not yet spread to the southern part. However, it is precisely this southern part that attracts our attention because it has been shown that seismic waves are amplified at a low frequency in this area. Within the framework of ongoing French–Ecuadorian research projects, a temporary network was installed for one year in order to carry out seismic interferometry based on ambient vibration measurements. The first results reveal a different structure in the north and south, which could explain the site effects obtained in the southern part

of the city (Pacheco *et al.* 2022). However, because the subsurface velocity model is not sufficiently known to compute numerical Green’s functions, as is the case for example for Mexico D.F. (Cruz-Atienza *et al.* 2016), and because we want to model the site effect reliably, we use an empirical Green’s function approach. We choose a stochastic summation method of empirical Green’s functions already proven in many cases and which requires few input parameters.

It has been shown that the stress drop (or more precisely the combined effect of stress drop and rupture velocity) is the parameter that mainly affects the acceleration peaks (Oth *et al.* 2017). We propose an original way to constrain the input stress drop using an equation giving the STF duration (mean and variability) as a function of the seismic moment from a global STF database, SCARDEC (Courboulex *et al.* 2016; Vallée & Douet 2016). These values are used as input for the empirical Green’s function summation method and make it possible to obtain a large number of synthetic seismograms representative of different source scenarios. This global approach is supported by the relative simplicity of the shapes of the source time functions at different magnitude scales, at least in the subduction zones (Meier *et al.* 2017).

The method we use is simple and quick to compute. The only part that needs to be undertaken with great care is the selection and analysis of small events used as EGFs. It is desirable to have
<figure><img src="image_11.png" /><figcaption>**Figure 9.** Simulation results ($M_w$ 7.8) for all the Quito stations and the four EGFs for PGA, PSA at period $T = 0.6$ s and $T = 2.8$ s. For comparison, the BC-hydro Abrahamson *et al.* (2016) GMMs' prediction is also presented (for Vs30 = 400 m s⁻¹). The x-axis represents the distance between each station and the rupture (Fig. 4).</figcaption></figure>

at least two of them (here we have four) to be sure not to choose a too particular small earthquake as EGF, which could bias the simulations. It is best to have earthquakes at different depths and different positions on the interface.

The results we obtained for the simulation of an earthquake of magnitude 7.8 showed that, without any adjustment of parameters, it was possible to obtain realistic synthetic seismograms. The comparison with the real records of the Pedernales earthquake of 2016 was a very good way for us to verify this. It showed us that, in this configuration (distance and situation of Quito regarding the subduction zone) the four small earthquakes we have used as EGFs gave very similar results, with no need for deeper or shallower events.

Similarly, the strong directivity of the Pedernales earthquake (southward, see Nocquet *et al.* 2016) had no visible influence on the Quito recordings and thus did not need to be considered.

Therefore, for this first paper on ground motion simulations in Quito, we did not find it useful to use more complex source models. We could have used a double corner frequency model as proposed by Archuleta & Ji (2016) and the code proposed by Nino *et al.* (2018), or the one generalized by Sedaghati *et al.* (2020) to describe the source. However, this would involve determining the values of the second corner frequency and there is currently no physical way to do so. It is clear that, for the $M_w$ 8.5 events, a refined extended source method that examines different scenarios with slip patches at different depths and with a non-constant rupture velocity, for example, is necessary in the future and will complement and specify the study proposed in this paper. Hybrid methods like the one presented by Castro-Cruz *et al.* (2021) could be an alternative.

The use of EGFs implies a linear hypothesis, that is the results obtained for small earthquakes should be similar to those obtained for large ones. The nonlinear effects of soils induced by strong vibrations have been studied in different contexts: Bonilla *et al.* (2011), Regnier *et al.* (2013) and Castro-Cruz *et al.* (2020) worked on the nonlinear behaviour of soils in Japan before, during and after the Tohoku, $M_w$ 9 earthquake. They found that high frequencies tended to decrease from a certain PGA level but low frequencies tended to be amplified. In the case of our study, the great earthquake being far from Quito, the level of the high frequencies is rather low (the maximum value of PGA is 0.2 g), so the nonlinear effects should also be weak for the high frequencies. On the contrary, for low frequencies like for Tohoku earthquake (Regnier *et al.* 2013), it is possible that the nonlinear effects will increase its values, resulting in even higher low-frequency values towards the south of the basin. The nonlinear effect due to a strong subduction earthquake would therefore be in the direction of an increase in the low-frequency levels.

Let us now examine the results obtained and try to determine whether an $M_w$ 8.5 event on the subduction zone constitutes a threat to the city of Quito. Concerning high frequencies (PGA values), the simulations we obtain are, most of the time, very close to those predicted by AB2016 (Fig 10): an average PGA value of 0.1 g and the +1 sigma value less than 0.2 g. This value is rather moderate but still could lead to significant damage to poorly constructed and thus vulnerable buildings, which are numerous in Quito. We compare the results at four stations with the design response spectra given by the National Ecuadorian Construction code (Norma Ecuatoriana de la Construcción, MIDUVI, 2015), and found that an $M_w$ 8.5 scenario should not drive any heavy damage to the most recent constructions (Fig. 12), if they are built respecting the code. Indeed, at 24MA, CIRC and LILI stations, located in the Quito basin, taking into account the possible soil class according to the NEC-15 assessed by Aguiar *et al.* (2016) of type C or D, corresponding to Vs30 ranging between 180 and 760 m s⁻¹, the design response spectra is much higher than the simulated one, on most of the entire period band. However, some of the simulations overpass the design spectra at station LILI (and the other stations of the south of Quito) at periods greater than 2 s. At ALLO station, set up on a site belonging to soil class A to B (Vs30 > 760 m s⁻¹), almost half of the simulations are overpassing the design spectra between 0.5 and 0.8 s. It is also noticeable that for this station at least, the maximum of the simulated response spectrum is found in a higher period band than the one corresponding to the plateau of the design spectrum.

Buildings of 7 to 12 stories [if we take into account the formula established by Perrault *et al.* (2020) for Quito's building stock,
<figure><img src="image_12.png" /><figcaption>**Figure 10.** Simulation results for $M_w$ 8.2 and 8.5 earthquakes at the Quito stations for PGA, PSA at period $T = 0.1$ s, $T = 0.6$ s and $T = 2.8$ s. Two scenarios are presented: one is centred on EGFS1 and the other one on EGF-N1. For comparison, the BC-hydro Abrahamson *et al.* (2016) GMMs' predictions are also presented (for Vs30 = 400 m s⁻¹).</figcaption></figure>
<figure><img src="image_13.png" /><figcaption>**Figure 11.** Examples of synthetic accelerograms ($M_w$ 8.5) obtained for the NS component of four stations in the Quito basin (24MA, ALLO, CIRC and LILI). The value to the right of each trace indicates the duration of the eSTF used for each simulation. The EGF-S1 derived from the earthquake of 2016 April 22 was used.</figcaption></figure>

$N = 14.3$ T], constructed on rocky outcrop in Quito, at least around ALLO site, could thus be at risk in this $M_w$ 8.5 scenario. Similarly, in the south of Quito, very high-rise building (buildings with more than 28 floors) could also potentially suffer from heavy damage.

Up to now, such buildings did not exist in the southern part of Quito. Nevertheless, in the future, the city will probably expand southwards and, given the pressure due to the constant increase in population, such high constructions could be proposed. Special care with the construction of this type of building is therefore highly recommended.

It is important to note that the predicted values of this study should be taken with prudence since the method we use does not properly take into account the complexity of very large earthquakes. Kurahashi & Irikura (2011) showed that during the great Tohoku ($M_w$ 9) earthquake, most of the high frequencies were generated at a depth greater than 35 km, in regions called SMGE (Strong Motion Generation Areas). Lay et al. (2012) report that it was the same for the Sumatra-Andaman ($M_w$ 9.2) and 2010 Chile ($M_w$ 8.8) megathrust earthquakes. Because the earthquakes we used as EGFs occurred at shallow depth between 10 and 15 km, we were not able to take into account this frequency-depth dependence. It is

clear that a complex simulation that takes into account this effect should be undertaken to well anticipate the high-frequency ground motions. It would also be crucial to take into account the fact that a huge earthquake could also break the surface and generate a large amount of surface waves, which would certainly enhance the low-frequency part of the signal.

## 7 CONCLUSION

We propose a simple approach where the input variability of the source time function is taken from a global database and used in a stochastic empirical Green's function method. With this method, we estimate the ground motion that could be generated in Quito by large subduction earthquakes. Our results show that the high frequencies generated by such earthquakes should be moderate but that, especially in the southern part of the city, we expect large amplitudes at low frequencies (around 0.3 Hz). We found that the building code of Quito is well adapted to resist to such earthquakes. It is therefore crucial that future constructions, especially tall buildings in the south of the city, strictly follow the imposed code in order not to suffer serious damage. However, the results obtained
<figure><img src="image_14.png" /><figcaption>**Figure 12.** Simulation results ($M_w$ 8.5) at four stations and two different EGFs represented as elastic response Spectra (damping 5 per cent) in accelerations (NS component). Mean values of the simulations are in a continuous line and mean value ± sigma are represented by dashed lines. GMMs' predictions (AB2016) are also presented (for Vs30 = 400 m s⁻¹). The design spectra of Quito are also plotted for four types of soil.</figcaption></figure>

with this simple approach need to be completed by a large number of extended fault simulations that properly take into account the complexity of large earthquakes in subduction zones.

## ACKNOWLEDGMENTS

We deeply thank the persons in charge of the seismic network and more generally all the researchers, engineers and students working at the IG-EPN of Quito. We thank Diego Mercerat and Daniel Pacheco for all the work carried out together on the Quito basin and fruitful discussions on the seismic hazard there. We thank Allie Hutchinson and Anne-Marie Duval for their readings and comments, Alain Dujardin and Laetitia Honoré for the previous work on the method. Thanks to Mario Ordaz, two anonymous reviewers and the associate editor Victor Cruz for helping us to strongly improve this manuscript.

This work was supported by the ANR-15-CE04-004 REMAKE project, the UCA/JEDI ANR-15-IDEX-01 project, the Geophysical Institute in Quito, the CEREMA and the Géoazur laboratory. It has been carried out within the framework of the Joint International Laboratory Séismes & Volcans dans les Andes du Nord (IRD LMI SVAN).

## DATA AVAILABILITY

The source time functions used to calibrate the input parameters can be found at http://scardec projects sismo ipgp fr/, last accessed: 2021 October 20.

The Python code used to build the aSTF can be found at https://gi tlab.com/dac/empirical-green-function, last access: 2021 October 20. Seismological data are available on request at https://www.igep n.edu.ec/eq20160416-data (Singaucho et al. 2016).

## REFERENCES

* Abrahamson, N., Gregor, N. & Addo, K., 2016. BC Hydro ground motion prediction equations for subduction earthquakes, *Earthq. Spectra*, **32**, 23–44.

* Aguiar, R., Rivas-Medina, A., Caiza, P. & Quizanga, D., 2016. Control spectra for Quito, *Nat. Hazards Earth Syst. Sci. Discuss.*, 1–19, doi:10.5194/nhess-2016-33

* Agurto-Detzel, H. et al., 2019. Ridge subduction and afterslip control after-shock distribution of the 2016 $M_w$ 7.8 Ecuador earthquake, *Earth planet. Sci. Lett.*, **520**, 63–76.

* Alvarado, A. et al., 2014. Active tectonics in Quito, Ecuador, assessed by geomorphological studies, GPS data, and crustal seismicity, *Tectonics*, **33**, 67–83.
Anderson, J., Bodin, P., Brune, J., Prince, J., Singh, S., Quaas, R. & Onate, M., 1986. Strong ground motion from the Michoacan, Mexico, earthquake, *Science*, **233**, 1043–1049.

Archuleta, R.J. & Ji, C., 2016. Moment rate scaling for earthquakes $3.3 \le M \le 5.3$ with implications for stress drop, *Geophys. Res. Lett.*, **43**, 12004–12011.

Beauval, C. et al., 2013. An earthquake catalog for seismic hazard assessment in Ecuador, *Bull. seism. Soc. Am.*, **103**, 773–786.

Beauval, C. et al., 2017. Comparison of observed ground-motion attenuation for the 16 April 2016 $M_w$ 7.8 Ecuador megathrust earthquake and its two largest aftershocks with existing ground-motion prediction equations, *Seismol. Res. Lett.*, **88**, 287–299.

Beauval, C. et al., 2018. A new seismic hazard model for Ecuador, *Bull. seism. Soc. Am.*, **108**, 1443–1464.

Bonilla, L.F., Tsuda, K., Pulido, N., Régnier, J. & Laurendeau, A., 2011. Nonlinear site response evidence of K-NET and KiK-net records from the 2011 off the Pacific coast of Tohoku earthquake, *Earth Planets Space*, **63**, 785–789.

Brune, J.N., 1970. Tectonic stress and the spectra of seismic shear waves from earthquakes, *J. geophys. Res.*, **75**, 4997–5009.

Castro-Cruz, D., Gatti, F. & Lopez-Caballero, F., 2021. High-fidelity broadband prediction of regional seismic response: a hybrid coupling of physics-based synthetic simulation and empirical Green functions, *Nat. Hazards*, **108**, 1997–2031.

Castro-Cruz, D., Régnier, J., Bertrand, E. & Courboulex, F., 2020. A new parameter to empirically describe and predict the non-linear seismic response of sites derived from the analysis of Kik-Net database, *Soil Dyn. Earthq. Eng.*, **128**, 105833.

Causse, M. & Song, S.G., 2015. Are stress drop and rupture velocity of earthquakes independent? Insight from observed ground motion variability, *Geophys. Res. Lett.*, **42**, 7383–7389.

Chalumeau, C. et al., 2021. Repeating earthquakes at the edge of the afterslip of the 2016 Ecuadorian $M_w$ 7.8 Pedernales earthquake, *J. geophys. Res.*, **126**, e2021JB021746.

Chatelain, J.-L. et al., 1999. Earthquake risk management pilot project in Quito, Ecuador, *GeoJournal*, **49**, 185–196.

Chounet, A. & Vallée, M., 2018. Global and interregion characterization of subduction interface earthquakes derived from source time functions properties, *J. geophys. Res.*, **123**, 5831–5852.

Cotton, F., Archuleta, R. & Causse, M., 2013. What is sigma of the stress drop?, *Seismol. Res. Lett.*, **84**, 42–48.

Courboulex, F., Converset, J., Balestra, J. & DelLouis, B., 2010. Ground-motion simulations of the 2004 $M_w$ 6.4 Les Saintes, Guadeloupe, earthquake using ten smaller events, *Bull. seism. Soc. Am.*, **100**, 116–130.

Courboulex, F., Vallée, M., Causse, M. & Chounet, A., 2016. Stress drop variability of shallow earthquakes extracted from a global database of source time functions, *Seismol. Res. Lett.*, **87**, 912–918.

Cruz-Atienza, V.M., Tago, J., Sanabria-Gómez, J.D., Chaljub, E., Etienne, V., Virieux, J. & Quintanar, L., 2016. Long duration of ground motion in the Paradigmatic Valley of Mexico, *Sci. Rep.*, **6**, doi:10.1038/srep38807.

Del Pino, I. & Yepes, H., 1990. *Apuntes para una historia sísmica de Quito*. Centro Histórico de Quito: Problemática y Perspectivas, Serie Quito. I. Municipio de Quito y Junta de Andalucía. Min. de Asuntos Exteriores de España, Quito, 67–100.

Douglas, J. & Aochi, H., 2008. A survey of techniques for predicting earthquake ground motions for engineering purposes, *Surv. Geophys.*, **29**, 187–220.

Font, Y., Segovia, M., Vaca, S. & Theunissen, T., 2013. Seismicity patterns along the Ecuadorian subduction zone: new constraints from earthquake location in a 3-D a priori velocity model, *Geophys. J. Int.*, **193**, 263–286.

Gallović, F., 2017. Azimuthal dependence of the ground motion variability from scenario modeling of the 2014 $M_w$ 6.0 South Napa, California, earthquake using an advanced kinematic source model, *Pure appl. Geophys.*, **174**, 3467–3478.

Godano, M., Bernard, P. & Dublanchet, P., 2015. Bayesian inversion of seismic spectral ratio for source scaling: Application to a persistent multiplet in the western Corinth rift, *Journal of Geophysical Research: Solid Earth*, **120**(11), 7683–7712.

Guéguen, P., Chatelain, J.-L., Guillier, B. & Yepes, H., 2000. An indication of the soil topmost layer response in Quito (Ecuador) using noise H/V spectral ratio, *Soil Dyn. Earthq. Eng.*, **19**, 127–133.

Gutenberg, B.U. & Richter, C.F., 1954. *Seismicity of the Earth and Associated Phenomena*, Read Books Ltd.

Hartzell, S. H., 1978. Earthquake aftershocks as Green's functions, *Geophysical Research Letters*, 5(), 1–4.

Honoré, L., Courboulex, F. & Souriau, A., 2011. Ground motion simulations of a major historical earthquake (1660) in the French Pyrenees using recent moderate size earthquakes: simulation of a major Pyrenean earthquake, *Geophys. J. Int.*, **187**, 1001–1018.

Joyner, W.B. & Boore, D.M., 1986. On simulating large earthquakes by Green's function addition of smaller earthquakes, *Earthq. Source Mech.*, **37**, 269–274.

Kanamori, H. & McNally, K.C., 1982. Variable rupture mode of the subduction zone along the Ecuador-Colombia coast, *Bull. seism. Soc. Am.*, **72**, 1241–1253.

Kelleher, J. A., 1972. Rupture zones of large South American earthquakes and some predictions, *Journal of Geophysical Research*, 77(), 2087–2103.

Kendrick, E., Bevis, M., Smalley, R., Brooks, B., Vargas, R.B., Lauría, E. & Fortes, L.P., 2003. The Nazca–South America Euler vector and its rate of change, *J. South Am. Earth Sci.*, **16**, 125–131.

Kohrs-Sansorny, C., Courboulex, F., Bour, M. & Deschamps, A., 2005. A two-stage method for ground-motion simulation using stochastic summation of small earthquakes, *Bull. seism. Soc. Am.*, **95**, 1387–1400.

Kurahashi, S. & Irikura, K., 2011. Source model for generating strong ground motions during the 2011 off the Pacific coast of Tohoku earthquake, *Earth Planets Space*, **63**, 571–576.

Laurendeau, A. et al., 2017. Low-frequency seismic amplification in the Quito Basin (Ecuador) revealed by accelerometric recordings of the RENAC network, *Bull. seism. Soc. Am.*, **107**, 2917–2926.

Lay, T. et al., 2012. Depth-varying rupture properties of subduction zone megathrust faults: megatrust rupture domains, *J. geophys. Res.*, **117**, doi:10.1029/2011JB009133.

León-Ríos, S. et al., 2019. 1D-velocity structure and seismotectonics of the Ecuadorian margin inferred from the 2016 Mw 7.8 Pedernales aftershock sequence, *Tectonophysics*, **767**, 228165.

Lorenzo, G.W.F., Santisi d'Avila, M.P.S., Deschamps, A., Bertrand, E., Mercerat, E.D., Foundotos, L. & Courboulex, F., 2018. Numerical and empirical simulation of linear elastic seismic response of a building: the case of Nice Prefecture, *Earthquake Spectra*, **34**, 169–196.

Madariaga, R., 1976. Dynamics of an expanding circular fault, *Bull. seism. Soc. Am.*, **66**, 639–666.

Mariniere, J. et al., 2020. Geodetic evidence for shallow creep along the Quito fault, Ecuador, *Geophys. J. Int.*, **220**, 2039–2055.

Meier, M.A., Ampuero, J.-P. & Heaton, T.H., 2017. The hidden simplicity of subduction megathrust earthquakes, *Science*, **357**, 1277–1281.

Meltzer, A. et al., 2019. The 2016 M w 7.8 Pedernales, Ecuador, earthquake: rapid response deployment, *Seismol. Res. Lett.*, **90**, 1346–1354.

Niño, M., Ayala, G. & Ordaz, M., 2018. Ground-motion simulation by the empirical Green's function method with a source defined by two corner frequencies and a two-stage summation scheme, *Bull. seism. Soc. Am.*, **108**, 901–912.

Nocquet, J-M. et al., 2014. Motion of continental slivers and creeping subduction in the northern Andes, *Nat. Geosci.*, **7**, 287.

Nocquet, J.-M. et al., 2017. Supercycle at the Ecuadorian subduction zone revealed after the 2016 Pedernales earthquake, *Nat. Geosci.*, **10**, 145–149.

Ordaz, M., Arboleda, J. & Singh, S.K., 1995. A scheme of random summation of an empirical Green's function to estimate ground motions from future large earthquakes, *Bull. seism. Soc. Am.*, **85**, 1635–1647.

Oth, A., Miyake, H. & Bindi, D., 2017. On the relation of earthquake stress drop and ground motion variability, *Journal of Geophysical Research: Solid Earth*, **122**(7), 5474–5492.

Pacheco, D., Mercerat, E.D., Courboulex, F., Bonilla, L.F., Laurendeau, A. & Alvarado, A., 2022. Profiling the Quito basin (Ecuador) using seismic ambient noise, *Geophys. J. Int.*, doi:10.1093/gji/ggab408.
Perrault, M., Guéguen, P., Parra, G. & Sarango, J., 2020. Modification of the data-driven period/height relationship for buildings located in seismic-prone regions such as Quito (Ecuador), *Bull. Earthq. Eng.*, **18**, 3545–3562.

Pulido, N., Yoshimoto, M. & Sarabia, A.M., 2020. Broadband wavelength slip model of the 1906 Ecuador–Colombia megathrust-earthquake based on seismic intensity and tsunami data, *Tectonophysics*, **774**, 228226.

Regnier, J., Cadet, H., Bonilla, L.F., Bertrand, E. & Semblat, J.-F., 2013. Assessing nonlinear behavior of soils in seismic site response: statistical analysis on KiK-net strong-motion data, *Bull. seism. Soc. Am.*, **103**, 1750–1770.

Reyes Once, M.F., 2020. Caracterización de la respuesta Sísmica de la Cuenca de Quito aplicando el Método de las componentes horizontales a vertical (H/V) a través del Ruido Ambiental, Quito.

Rodgers, A.J., Pitarka, A., Petersson, N.A., Sjögreen, B. & McCallen, D.B., 2018. Broadband (0-4 Hz) ground motions for a magnitude 7.0 Hayward Fault earthquake with three-dimensional structure and topography: an M 7 Hayward Fault earthquake, *Geophys. Res. Lett.*, **45**, 739–747.

Rolandone, F. et al., 2018. Areas prone to slow slip events impede earthquake rupture propagation and promote afterslip, *Sci. Adv.*, **4**, eaao6596.

Salichon, J., Kohrs-Sansorny, C., Bertrand, E. & Courboulex, F., 2010. A M<sub>w</sub> 6.3 earthquake scenario in the city of Nice (southeast France): ground motion simulations, *J. Seismol.*, **14**, 523–541.

Sedaghati, F., Tavakoli, B. & Pezeshk, S., 2020. A generalization of the stochastic summation scheme of small earthquakes to simulate strong ground motions, *Pure appl. Geophys.*, **177**, 3713–3732.

Singaucho, J. C., Laurendeau, A., Viracucha, C. & Ruiz, M., 2016. Observaciones del sismo del 16 de Abril de 2016 de magnitud M<sub>w</sub> 7.8, Intensidades y Aceleraciones. From http://www.igepn.edu.edu/1324-info rme-sismico-especial-n-18-2016.html (accessed: 20 April 2021).

Tsang, L. et al., 2019. Imaging rapid early afterslip of the 2016 Pedernales earthquake, *Earth planet. Sci. Lett.*, **524**, 115724.

Ulrich, T., Gabriel, A.A., Ampuero, J.-P. & Xu, W., 2019. Dynamic viability of the 2016 M<sub>w</sub> 7.8 Kaikōura earthquake cascade on weak crustal faults, *Nat. Commun.*, **10**, doi:10.1038/s41467-019-09125-w

Vaca, S., Vallée, M., Nocquet, J.-M. & Alvarado, A., 2019. Active deformation in Ecuador enlightened by a new waveform-based catalog of earthquake focal mechanisms, *J. South Am. Earth Sci.*, **93**, 449–461.

Vallée, M., 2013. Source time function properties indicate a strain drop independent of earthquake depth and magnitude, *Nat. Commun.*, **4**, doi:10.1038/ncomms3606

Vallée, M., Charléty, J., Ferreira, A.M.G., Delouis, B. & Vergoz, J., 2011. SCARDEC: a new technique for the rapid determination of seismic moment magnitude, focal mechanism and source time functions for large earthquakes using body-wave deconvolution: wave deconvolution and earthquake parameters, *Geophys. J. Int.*, **184**, 338–358.

Vallée, M. & Douet, V., 2016. A new database of source time functions (STFs) extracted from the SCARDEC method, *Phys. Earth planet. Inter.*, **257**, 149–157.

Vieux-Champagne, F. et al., 2017. Experimental analysis of a shake table test of timber-framed structures with stone and earth infill, *Earthq. Spectra*, **33**, 1075–1100.

Wang, H., Ren, Y., Wen, R. & Xu, P., 2019. Investigating the contribution of stress drop to ground-motion variability by simulations using the stochastic empirical Green’s function method, *Pure appl. Geophys.*, **176**, 4415–4430.

Wang, H., Wen, R. & Ren, Y., 2017. Simulating ground-motion directivity using stochastic empirical Green’s function method, *Bull. seism. Soc. Am.*, **107**, 359–371.

Wennerberg, L., 1990. Stochastic summation of empirical Green’s functions, *Bull. seism. Soc. Am.*, **80**, 1418–1432.

Yamanaka, Y., Tanioka, Y. & Shiina, T., 2017. A long source area of the 1906 Colombia–Ecuador earthquake estimated from observed tsunami waveforms, *Earth Planets Space*, **69**, doi:10.1186/s40623-017-0750-z

Ye, L., Kanamori, H., Avouac, J.-P., Li, L., Cheung, K.F. & Lay, T., 2016. The 16 April 2016, M<sub>w</sub> 7.8 (M<sub>s</sub> 7.5) Ecuador earthquake: a quasi-repeat of the 1942 M<sub>s</sub> 7.5 earthquake and partial re-rupture of the 1906 M<sub>s</sub> 8.6 Colombia Ecuador earthquake, *Earth planet. Sci. Lett.*, **454**, 248–258.

Yepes, H., Audin, L., Alvarado, A., Beauval, C., Aguilar, J., Font, Y. & Coton, F., 2016. A new view for the geodynamics of Ecuador: implication in seismogenic source definition and seismic hazard assessment, *Tectonics*, 2015TC003941, doi:10.1002/2015TC003941

Yoshimoto, M. et al., 2017. Depth-dependent rupture mode along the Ecuador–Colombia subduction zone: depth dependent rupture mode, *Geophys. Res. Lett.*, **44**, 2203–2210.

## DETAILS ON THE PROCEDURE

A seismogram at any point can be expressed as

$$S(f) = \text{STF}(f) \times \text{Path}(f) \times \text{Site}(f). \quad (A1)$$

If the system is linear, for a strong and a weak earthquake that comes from the same fault (same distance and same focal mechanism), the path effects and site effects are the same. Then the scaling function R can be written as

$$R(f) = \frac{S_{\text{strong}}}{S_{\text{weak}}} = \frac{\text{STF}_{\text{strong}}}{\text{STF}_{\text{weak}}}. \quad (A2)$$

If the two earthquakes follow the ω⁻²-model (Brune 1970):

$$R(f) = \frac{M_o}{m_o} \cdot \frac{1 + \left(\frac{f}{F_c}\right)^2}{1 + \left(\frac{f}{F_c}\right)^2}, \quad (A3)$$

where M<sub>o</sub> and F<sub>c</sub> are the seismic moment and the corner frequency of the strong earthquakes (the one to be simulated, called the target event) and m<sub>o</sub> and f<sub>c</sub> are the seismic moment and the corner frequency of the weak event taken as empirical Green’s function (and called EGF).

With the previous conditions, Fig. A1 shows the general shape of the scale between two earthquakes (R). At low frequencies, R depends exclusively on the ratio of the seismic moment between both earthquakes. At higher frequencies, the stress drop ratio between both earthquakes (represented by the ratio of the corner frequencies f<sub>c</sub> and F<sub>c</sub> in this figure) also influences the shape of R.

<figure><img src="image_15.png" /><figcaption>Figure A1. General shape of the scale function R. The axes are in log₁₀ scale. F<sub>c</sub> and f<sub>c</sub> represent respectively the corner frequencies of the strong target event and the weak event taken as EGF. M<sub>o</sub>/m<sub>o</sub> is the ratio between the seismic moments of the strong and weak events.</figcaption></figure>

Eq. (A3) represents the ratio of a theoretical case. However, to create a realistic result, the variability of the release of energy in the source (for a constant value of M<sub>o</sub>) must be considered. Joyner & Boore (1986) proposed to sum lags of time to represent the spatial and random release of energy during the source. These random
delays are generated in a way that in the frequency domain eq. (A3) is on average satisfied.

Several methods were proposed to generate these delays. Wennerberg (1990) proposed a method where the time delays are generated by groups in a two-stage method. This way of generation makes it possible to obtain realistic and different histories that can be associated with the real ones. This method avoids the convergence when a large number of time delays are done in the simulations. For this reason, this EGF method is applicable in simulations of stronger events, obtaining a large range of different possible rupture processes. With this model of two-step generation, the ratio R is generated with a previous step in the time domain. A time function (r) is defined as

$$r^{(i)} = \kappa \sum_{d=0}^{\eta_d} \left[ \sum_{c=0}^{\eta_c} \delta(t - t_c - t_d) \right], \quad (A4)$$

where $r^{(i)}$ is composed of scaled points that are also composed of a summation of impulses that are distributed in time randomly. The index i represents the index of the simulation, and the factor $\kappa$ is a coefficient of scaling. $\delta$ refers to Dirac's delta function, where $t_c$ and $t_d$ are random delays in time. $\eta_d$ and $\eta_c$ are the numbers of pulses in which the signal is divided.

The Dirac's function ($\delta$) in eq. (A4) cannot be defined, and a discrete approximation is necessary since all the process is done with numerical methods. The function $\delta$ is approximated as

$$\delta(t - t_c - t_d) \approx \frac{1}{\Delta t} \cdot \hat{\delta} \left( \text{int} \left( \frac{t - t_c - t_d}{\Delta t} \right), 0 \right). \quad (A5)$$

The term $\Delta t$ represents the time step used to store the accelerogram. This term is necessary because the discrete approximation is done with the Kronecker delta ($\hat{\delta}$), so to guarantee that the energy of $r^{(i)}$ is not affected, this term must be added. The function int is a truncation process to remove the decimals of the number in the interior.

The method of EGF that is used (Wennerberg 1990; Kohrs-Sansorny et al. 2005) makes the equality in eq. (A3) to be accomplished on average for all the simulations. To achieve this, the scale ($\kappa$), the number of delays ($\eta_c$ and $\eta_d$), and the distributions of probability for $t_c$ and $t_d$ are fitted. $\eta_d$ and $\eta_c$ are defined as:

$$\eta = \eta_d \cdot \eta_c \quad \text{where: } \eta = \left( \frac{f_c}{F_c} \right)^4 \quad (A6)$$

Eq. (A6) guarantees that there are enough impulses to estimate the signal at high frequencies, but also that there are not too many pulses to make that the simulations tend to be unrealistic and identical between them. The coefficient $\kappa$ is obtained to scale the number of used points ($\eta$) with the scale function that is searched (eq. A3).

$$\kappa = \frac{M_o}{m_o \cdot \eta} \quad (A7)$$

The delay times $t_c$ and $t_d$ (eq. A5) are distributed with a density probability function that makes, in average, the Fourier transform amplitude of $r^{(i)}$ will be equal to the expected scale R (eq. A3).

To compute this probability density function for the time lags, we follow the method of Kohrs-Sansorny et al. (2005) that combines the density of probability of Ordaz et al. (1995) and the two-step summation of Wennerberg (1990). The probability density function for the time delays is

$$\rho = \int_{-\infty}^{\infty} \frac{\sqrt{1 + \left(\frac{2}{1+(F_2/F_1)^2}\right) \cdot (f/F_1)^2}}{1 + (f/F_2)^2} \cdot e^{i2\pi ft} df, \quad (A8)$$

where for $\rho_c$, the density function of $t_c$, the values $F_1$ and $F_2$ are $F_d$ and $F_c$, respectively. For $\rho_d$, the density function of $t_d$, the values $F_1$ and $F_2$ are $F_d$ and $F_c$. $F_c$ and $f_c$ being the corner frequencies of the strong earthquake and the weak earthquake respectively. $F_d$ is defined as

$$F_d = \eta_c^{1/4} \cdot F_c. \quad (A9)$$

Generating the time delays $t_c$ and $t_d$ with adequate distributions defined in eqs (A6) and (A9) allows us to find the ratio between the spectra from weak and strong ground motions (eq. A2)

$$S_{\text{strong}}(f) = S_{\text{weak}}(f) \cdot \hat{R}^{(i)}(f), \quad (A10)$$

where $S_{\text{strong}}$ is the amplitude spectrum (FAS) of the signal at the surface of the target signal. $S_{\text{weak}}$ the Fourier spectrum of the weak record (EGF). $\hat{R}^{(i)}$ is the Fourier transform of the summation of the impulses defined in eq. (A4). Eq. (A10) in time domain represents a convolution defined as:

$$s_{\text{strong}}^{(i)}(t) = \text{EGF}(t) * r_i(t) \quad (A11)$$

where $s_{\text{strong}}^{(i)}$ is one of the possibilities of the accelerogram of the target earthquake, EGF is the accelerogram used from the weak ground motion, and $r_i$ is the scale founded with the eq. (A4).

The scaling function $R(f)$ (eq. A3), at low frequencies depends only on the seismic moment ratio. It means that the predicted earthquake at low frequencies is not affected by the variation of the corner frequency of the target earthquake. This independence of the target event to the corner frequency occurs before $f$ approximates to 0.5 $F_c$.

At high frequencies, $R(t)$ also trends to a fixed value:

$$R = \frac{M_o}{m_o} \cdot \left( \frac{F_c}{(f_c)} \right)^2 .$$

This means that high frequencies of the predicted earthquake will be affected by the corner frequency of the target event $F_c$. The definition of high frequencies is marked approximately when $f$ is higher than 2 $f_c$. In the method we propose, we directly take the variability of the stress drop through the duration of the STF (T). The corner frequency of the target event is defined as $F_c = q/T$, with $q$ as a constant value that depends on the geometry of the fault (see the annex of Godano et al. 2015). In this paper, $q$ is taken equal to 0.8.