the Creative Commons Attribution 4.0 License.
the Creative Commons Attribution 4.0 License.
Estimating ground motion intensities using simulationbased estimates of local crustal seismic response
Himanshu Agrawal
John McCloskey
It is estimated that 2 billion people will move to cities in the next 30 years, many of which possess high seismic risk, underscoring the importance of reliable hazard assessments. Current ground motion models for these assessments typically rely on an extensive catalogue of events to derive empirical ground motion prediction equations (GMPEs), which are often unavailable in developing countries. Considering the challenge, we choose an alternative method utilizing physicsbased (PB) ground motion simulations and develop a simplified decomposition of ground motion estimation by considering regional attenuation (Δ) and local site amplification (A), thereby exploring how much of the observed variability can be explained solely by wave propagation effects. We deterministically evaluate these parameters in a virtual city named Tomorrowville, located in a 3Dlayered crustal velocity model containing sedimentary basins, using randomly oriented extended sources. Using these physicsbased empirical parameters (Δ and A), we evaluate the intensities, particularly peak ground acceleration (PGA), of hypothetical future earthquakes. The results suggest that the estimation of PGA using the deterministic Δ−A decomposition exhibits a robust spatial correlation with the PGA obtained from simulations within Tomorrowville. This method exposes an orderofmagnitude spatial variability in PGA within Tomorrowville, primarily associated with the nearsurface geology and largely independent of the seismic source. In conclusion, advances in PB simulations and improved crustal structure determination offer the potential to overcome the limitations of earthquake data availability to some extent, enabling prompt evaluation of ground motion intensities.
 Article
(9400 KB)  Fulltext XML

Supplement
(25965 KB)  BibTeX
 EndNote

In the Global South, the absence of seismic catalogues impedes ground motion predictions that are crucial for earthquakeaware urban planning.

Physicsbased simulations can use hypothetical earthquakes to estimate ground motions without extensive earthquake data availability.

The primary source of shortscale variability in ground motion is the local subsurface geology, making it a crucial focal point.
The United Nations Human Settlements Programme (UNHabitat) forecasts that by 2050 some 2 billion new citizens will move to urban centres so that, by then, some 68 % of the world's population will live in cities (United Nations, 2022). It is estimated that 95 % of this urbanization will happen in the Global South. Urban population growth is often accommodated by rapid urban expansion in areas with welldocumented seismic risk. The problems of understanding and reducing disaster risk in such rapid development are significant, and while this expansion presents a major global challenge, it also provides a timelimited opportunity to provide evidencebased decision support for this new development (UNISDR, 2015). Efforts in earthquake risk reduction through urban planning guided by highresolution hazard assessment could reduce disaster risk for hundreds of millions of these future citizens. This approach also provides a costefficient method by concentrating on new constructions, where the expenses related to implementing effective earthquakeresistant design and construction are significantly lower compared to the costs of retrofitting at a later stage.
Seismic hazard analysis informs building codes, constraining the construction of new development in earthquakeprone areas through the development of ground motion models (Baker et al., 2021; Bradley, 2019; McGuire, 2008; Stirling et al., 2012; Stirling, 2014; Kramer and Mitchell, 2006; Kramer, 1996). Observed ground shaking is a result of the interaction between a range of individually heterogeneous fields and processes, leading to great complexity in even the simplest relationships. Measures of groundshaking intensity, for example, show an expected systematic decrease with the distance between the observation and source, but the systematics are overprinted by the interactions between the complexities of the event and the crustal volume explored by the seismicwave train. The result is highamplitude variability in the observed intensity. Note that the uncertainty in the observations, in either intensity or distance, makes only a small contribution to this variability; the variability is an intrinsic part of the process.
Consider a series of events recorded at a large number of sensors. In the commonly applied approach, the analyst chooses a functional form for the systematic decay of intensity and uses some fitting procedures to estimate its parameters. The resulting model is commonly known as a ground motion model (GMM) (Douglas and Edwards, 2016; Douglas and Aochi, 2008) and takes the form
where IM is the required intensity measure, μ_{lnIM} is the estimated meanfield intensity, σ_{lnIM} is an estimate of the variability around the mean that is usually assumed to conform to a lognormal distribution and ϵ is the standard normal variate.
It is important to note that the μ_{lnIM} term does not just describe the attenuation of intensity with distance. Common forms of μ_{lnIM} attempt to parameterize descriptions of the physics of the entire process, including source properties, such as the focal mechanism and the resulting directivity, as well as the local response of the site using estimates of V_{s30} (timeaveraged shearwave velocity in the top 30 m) and κ (highfrequency attenuation parameter), for example (Shi and Asimaki, 2017; Aki, 1993; Bradley, 2011; Kaklamanos et al., 2013; Hough and Anderson, 1988; Borcherdt and Glassmoyer, 1992). Expressions for μ_{lnIM} in current GMMs include numerous parameters, use advanced statistical techniques to fit these complex functions and represent a practical approach to a fundamentally intractable problem (Douglas and Edwards, 2016).
In practice, an ergodic assumption is invoked in GMM development by aggregating the data from multiple spatial locations that are assumed to be equivalent to the distribution in time (Anderson and Brune, 1999). However, with the increasing data for a particular tectonic area, the nonergodic or partial nonergodic approaches are favoured, modifying μ_{lnIM} and σ_{lnIM} based on calibration with the local data that are available (Stewart et al., 2017; Bradley, 2015; RodriguezMarek et al., 2014). It is observed that a major component of ground motion amplification can be associated with the local geological factors, e.g. sedimentary basins (Graves et al., 1998; Zhu et al., 2018; Pilz et al., 2011), surface topography (Maufroy et al., 2012; Wang et al., 2018; Lee et al., 2009) and soil conditions (Bazzurro and Cornell, 2004; Cramer, 2003; Torre et al., 2020). Hence, the general practice in GMM development is dominated using nearsurface sitespecific parameters (for example V_{s30} and κ). It is suggested that these nearsurface parameters might exhibit strong correlations with geological features at greater depths, like basin depth parameters (Z_{xx}) (Chiou and Youngs, 2014; Tsai et al., 2021; Kamai et al., 2016), and consequently the amplification. However, opposing studies show that the amplification patterns might not necessarily correlate with these parameters (Castellaro et al., 2008; Mucciarelli and Gallipoli, 2006; Pitilakis et al., 2019), for example, sites with velocity profiles which are not monotonically increasing with depth. This highlights the necessity to investigate more regional geological structure to better understand the complexities of ground motion amplification.
Recently, the advances in computational capabilities and understanding the physical processes have made it possible to use physicsbased (PB) simulations for modelling ground motions (Smerzini and Villani, 2012; Taborda and Roten, 2014; Bradley, 2019; Graves and Pitarka, 2010). PB simulations are carried out by numerical modelling of the entire process of rupture characterization and seismicwave propagation through the earth's potentially complex crust. However, the high computational cost and complex input requirements associated with them restrict the largescale usage of these methods, particularly in 3D. As a consequence the relative contribution of these processes to the total observed variability has been relatively unexplored compared to that of local shallow (decametre) site conditions.
Two immediate problems emerge in enacting the current approaches to ground motion modelling in the context of rapid urbanization in the Global South described above. Firstly, understanding ground motion requires extensive seismic databases recording appropriate measures of intensity from a large number of earthquakes recorded with a network of sensors in the area of interest, for example, the Pacific Earthquake Engineering Research Center Next Generation Attenuation (PEER NGA) databases (Ancheta et al., 2014; Spudich et al., 2013; Atkinson and Boore, 2006). Such catalogues necessitate the deployment of seismometers for many years, even in the most seismically active areas for which is not possible to address the current timecritical problem (Freddi et al., 2021). Secondly, urban development projects require hazard information at unusually high resolution. Urban flood modelling and landslide susceptibility estimates, for example, typically strive to use digital terrain models with 2 m resolution supplemented by highresolution geotechnical assessments (Jenkins et al., 2023). Seismic intensity also varies significantly over the scale of interest for urban planning, particularly where development is planned over sedimentary basins or near to coasts or rivers with strong spatial contrasts in subsurface seismic velocity (Bielak et al., 1999; see also, Cadet et al., 2011; Foti et al., 2019). Some efforts have been made to incorporate these factors into ground motion prediction equations (GMPEs; Chiou and Youngs, 2014; Campbell and Bozorgnia, 2014; Abrahamson et al., 2014; Marafi et al., 2017); however, the extensive information required to accurately characterize such effects remains a challenge. As a result, the potential for risk reduction with a high cost–benefit ratio that would accrue from a highresolution understanding of ground motion variability remains elusive. Typically, GMMs developed in datarich countries of the Global North are reconditioned for deployment in areas for which they have no obvious physical validity (Nath and Thingbaijam, 2011; Hough et al., 2016). At best, this leads to poor spatial resolution precluding the detailed site classification that is critical for seismic microzonation studies needed for costeffective urban planning (Ansal et al., 2010). The development of appropriate techniques for rapid, local, highresolution seismic hazard assessment is a significant global challenge.
In this research, we approach this challenge using a simplified decomposition of ground motions into parametric relations explaining the regional and local variations in the measured intensity. We focus on the effects only due to the sedimentary basins, which are known to enhance the amplitude and duration of seismic waves through frequencydependent focusing, trapping and resonance (Frankel, 1993; Yomogida and Etgen, 1993; Castellaro and Musinu, 2023). We demonstrate the usefulness of PB simulations in capturing the primary lowfrequency (LF), < 1 Hz, sedimentary basin effects that contribute to the variation in ground motion within an urban area situated within a seismically active region. We show, to first order, seismic intensity decays along the wave path according to the integrated rheological properties of the region, and this is concurrently subject to relative amplification specific to any point on the surface. We first provide the theoretical physical basis for the decomposition and then describe the simulation domain and the numerical scheme used to explore it. We then describe how the main elements of the problem, i.e. regional meanfield attenuation (Δ) and local site amplification (A) (explained in the subsequent section), can be extracted from the simulations and demonstrate their use in the reconstruction of originally simulated intensities. We highlight that the assessment of these reconstructed intensities is not notably influenced by source characteristics (such as location and directivity). Therefore, calibrating these parameters and understanding shortscale ground motion amplification variability can address the challenge posed by the lack of earthquake data. We suggest that this approach, when extended to include high frequency (HF), might provide an improved relative seismic risk assessment in the form of more reliable microzonation maps at the scale of urban planning, which is based on rapid seismological site characterization in the absence of longduration seismic catalogues.
Using the seismic representation theorem (De Hoop, 1958; Knopoff, 1956), in polar coordinates the displacement U_{δ,ε} recorded at site ε for pointsource earthquake δ is given by
where r is the distance between the source and receiver, θ and ∅ are the positional angles in a spherical coordinate system, f_{δ} is a force vector at δ, and G is the elastodynamic Green's function providing the displacement at ε due to f_{δ}. Since we consider the peak displacement in an elastic medium in what follows, this equation is time invariant.
Consider a receiver at point ε that experiences displacements due to sources of a given seismic moment at point δ (see Fig. 1). The average logarithm of the peak displacement field for all possible point sources δ_{r} at distance r from receiver ε can then be expressed as
where $\stackrel{\mathrm{\u203e}}{\mathrm{ln}\left({U}_{{\mathit{\delta}}_{r}\mathit{\epsilon}}\right)}$ then represents the expectation value for the intensity at ε due to all possible events at distance r. In this formulation, we consider point sources without any particular focal mechanism, so Eq. (3) might be considered an integration over all possible focal mechanisms at all possible points on the hemisphere.
Integrating over all receivers Ω_{ε} on the surface of the domain,
which then provides a meanfield estimate of the expected intensity for any source–receiver pair separated by distance r, and a graph of $\stackrel{\mathrm{\u203e}}{\mathrm{ln}\left({U}_{(\mathit{\delta}\mathit{\epsilon}{)}_{r}}\right)}$ against r represents the meanfield decay of intensity with distance throughout the entire volume.
The response at a particular location on the surface to any specific event at some distance r will, of course, be subject to the source, path and site effects, all contributing to some local modification of the meanfield expectation. Consider the ground motion at receiver ε due to any source δ; again, peak displacement (U_{δ,ε}) can be calculated using the representation theorem, this time giving
This peak ground displacement U_{δ,ε} varies with ε, but from Eq. (4), we know its mean across the surface is $\stackrel{\mathrm{\u203e}}{\mathrm{ln}({U}_{\left(\mathit{\delta}\mathit{\epsilon}{)}_{r}\right)}}$. Normalizing the U_{δ,ε} by $\stackrel{\mathrm{\u203e}}{\mathrm{ln}\left({U}_{(\mathit{\delta}\mathit{\epsilon}{)}_{r}}\right)}$ removes the meanfield decay leading to normalized displacement ${\widehat{U}}_{\mathit{\delta},\mathit{\epsilon}}$ given by
Finally, to encapsulate the effect of all possible sources at each receiver, this normalized displacement can be integrated for the entire source field (Ω_{δ}), giving
This $\stackrel{\mathrm{\u203e}}{\mathrm{ln}\left({\widehat{U}}_{\mathit{\epsilon}}\right)}$ describes a local normalized amplification expected at any point for all possible sources. This can be considered the integrated effect of the whole wave path from all possible sources that is dominated near ε, where these paths converge. This term introduces the empirical sitespecific variability using the normalized intensity of a suite of earthquakes of any magnitude.
Equations (4) and (7) now allow us to express the final estimate of the intensity measure as
For the sake of simplicity, for an event at i, observed at location j, separated by distance r, lnΔ_{r} is used to denote the first term, mean intensity decay $\stackrel{\mathrm{\u203e}}{\mathrm{ln}\left({U}_{{\left(\mathit{\delta}\mathit{\epsilon}\right)}_{r}}\right)}$, and lnA_{j} defines the second term describing amplification, $\stackrel{\mathrm{\u203e}}{\mathrm{ln}\left({\widehat{U}}_{\mathit{\epsilon}}\right)}$. Now, Eq. (8) can then be rewritten as
where IM_{ij} is a nonspecific intensity measure recognizing that the argument so far may be generalized to peak velocity or acceleration. IM_{ij} then provides an estimate of the intensity of ground motion based on the meanfield expected intensity at distance Δ_{r}, integrated over the entire crustal volume under consideration, and relative amplification A_{j} due to the integrated effect of the seismic velocity structure around the site. Both terms on the righthand side are properties of the crust, regional and local, and do not include extended descriptions of the earthquake source, as we show in the next section. Equation (9) defines the Δ−A decomposition, a static ground motion model that emphasizes local geology rather than the descriptions of the earthquake source.
In practice, mean field Δ and amplification A, can both be calibrated through simulationbased estimates for a given domain; hence the basis is essentially nonergodic, but it is different than databased statistically estimated parameters used in typical nonergodic GMMs (e.g. Landwehr et al., 2016; Kuehn et al., 2019). The spatial coefficients estimated in these nonergodic model are datadependent; hence in order to find potential drivers of GM variability in datasparse regions, there is very little scope to use these models. To clarify, the motivation for the potential use of the Δ−A method is to target the datasparse regions without the extensive availability of earthquake catalogues.
To explore the behaviour and stability of Δ and A (in Eq. 9) and how they might be estimated in practice, we use a virtual world that allows for the exploration of the ideas in the absence of uncertainty but which allows for the introduction of precisely constrained variability. We use a virtual crustal environment, as shown in Fig. 2a and b, that incorporates a simplified subsurface velocity structure centred on a shallow and a deep river basin overlying a crystalline basement to which simplified velocities have been assigned. The description of the domain includes depthvarying density (ρ), shear wave speed (V_{s}), primary wave speed (V_{p}) and anelastic attenuation factors (Q_{p},Q_{s}), and it is determined based on the assumed values of these parameters at the surface of the shallow basin (river channel), deep basin and basement (Brocher, 2005, 2008). The reader is referred to Sect. 3.1 of Jenkins et al. (2023) for detailed description of the crustal domain and earthquake moment distribution. Alternatively, this information is also accessible in the Supplement (Table S1 and Fig. S1).
In the middle of the crustal domain, we locate a virtual urban environment Tomorrowville (Menteşe et al., 2023; Gentile et al., 2022; Cremen et al., 2023; Wang et al., 2023; Jenkins et al., 2023). The geology of Tomorrowville is based on a stretch of the Nakhu River valley on the outskirts of Lalitpur to the south of Kathmandu though the velocity structure described here that extends far to the north and south and does not represent the actual subsurface seismic velocity in the area. Instead, we simply generate a hypothetical nearsurface velocity structure representative of any urban settlement located around a river channel set in a deeper and wider sedimentary basin. The depths of shallow and deep basins in Tomorrowville are presented in Fig. 2c and d.
The random distribution of 40 thrustfaulting earthquakes (EQ1 to EQ20 are M_{w} 6 and EQ21 to EQ40 are M_{w} 5) is simulated across the domain (see Fig. 2e, f) using an established physicsbased solver, SPEED, which uses spectral element method (SEM) for solving the wave propagation equations (Mazzieri et al., 2013; Paolucci et al., 2014; Smerzini et al., 2011). SEM combines the geometrical flexibility of the finiteelement method (FEM), i.e. the capability to naturally account for irregular interfaces and mesh adaptivity, with the high spectral accuracy, i.e. the exponential convergence rate to the exact solution that results in a fewer number of grid points per wavelength to maintain low dispersion. The crustal domain has a minimum shear wave velocity of 250 m s^{−1} and the smallest element size of 200 m with the spectral degree of 4; hence, the simulations are able to resolve for the vibrational periods greater than 0.8 s. Fault plane dimensions are determined using widely used empirical relationships developed by Wells and Coppersmith (1994). Kinematic characterization of rupture model is done based on the model developed by Liu et al. (2006) and Schmedes et al. (2013), in which the correlation between the slip, rise time, peak time and rupture velocity among the subfaults are derived based on a large ensemble of dynamic rupture simulations of dipping faults. The moment distribution remains same for each magnitude ensemble, but the strike and dip are varied. This distribution of rupture scenarios produces a wide range of expected source directivity for any location. The peak ground acceleration (PGA) maps shown in Fig. S2 and Movie S1 are referred to for the visualization of the source orientation and their corresponding effects across the surface of entire domain. The wavefront evolution for EQ1 can also be found in Movies S2, S3 and S4 of the Supplement as well.
The Δ−A decomposition, developed theoretically above (Sect. 2), includes no source variability, whereas any attempt to understand seismic hazard must. The azimuth of the events from the seismometer with respect to the dominant velocity anisotropy introduced by the river basin will also contribute to the expected ground motion variability. The aim of this paper is not to examine the influence of these features on the observed local intensity; that will follow in a later work. Instead, we simply explore the extent to which the relative amplification term, A_{j}, might act as a usable proxy that, to first order, governs the intensity variation across an urban area, irrespective of the source orientation. This might be considered a lower bound on the skill of Eq. (9) in providing the basis for a static sitedependent ground motion model that might be improved later by the introduction of a source term to be constrained by the structural fabric and stress state around any specific location.
The simulation results are used to estimate Δ for the crustal domain and A for Tomorrowville (Eq. 9). The geometric mean of horizontal components of PGA values are used as an intensity measure for all of the rupture scenarios.
To calculate Δ, we uniformly sample the surface of the crustal domain, which is a practical and computationally inexpensive approach to approximate the integration in Eq. (4). In the entire simulation domain, a random set of 100 recording locations is chosen (see green triangles in Fig. 3a), for which estimates of the PGA are simulated for every event, generating a large number of estimates of the peak amplitude for different epicentral distances giving the data points for events of magnitude 5 and 6 shown in Fig. 3b. We use simple leastsquares regression to the decay equation
where Δ_{r} is an estimation of the meanfield intensity measure Δ_{r} (introduced in Eq. 9); r is the epicentral distance; and a, b and c are the empirical parameters evaluated from the datafitting procedure which might be modified without loss of insight (Fig. 3b). The choice of 100 recording locations for Δ_{r} estimation can have inherent uncertainties based on the selection. For instance, if the stations are predominantly concentrated in the basin, it could result in higher intensities in Fig. 3b, consequently causing an upward shift in the meanfield curve. However, such a scenario would not uniformly sample the entire domain as intended; hence, the current choice of stations seems satisfactory.
It should be noted that the regression method chosen here does not distinguish the repeatable (withinevent) and nonrepeatable (betweenevents) effects, which is determined the fact that each source used here is characteristically similar and is recorded at the exact same set of receivers. Assuming the entire domain has a homogeneous earthquake distribution, each recording is considered independent, irrespective of whether the seismic energy originated from the same or different sources. The concept of earthquake source homogeneity implies that in a scenario with limited prior knowledge of the tectonics in the area, a reverse faulting earthquake could potentially occur at any azimuth with respect to the city.
We now must turn our attention to the variability in the data around the curves (Fig. 3b) and will focus on the Tomorrowville subdomain. Note that any numerical uncertainties due to the calculation, conditional on the input geological structure, are negligible compared to the variability observed in Fig. 3b. Hence, given the assumption that the simulation is providing accurate estimates in a virtual setting, each point in Fig. 3b accurately represents the local peak amplitude of waves from a particular event recorded at a single station. To estimate A_{j} for any location j, the PGA values from all events are extracted for the Tomorrowville domain (Fig. 4a). The linear interpolation of intensities is used to provide these highresolution maps, which sample Tomorrowville at an approximate grid spacing of 28 m.
As an example, the PGA from earthquake 1 (EQ1) is shown along with the spectral accelerations (5 % damped) at 10 stations, S1 to S10 (Fig. 4b, c). Please note that these receivers are positioned within the Tomorrowville domain and are not accounted for in the wider receiver distribution illustrated in Fig. 3a for the evaluation of Δ_{r}. It can be clearly seen that the basin area is showing strong amplification resulting in higher PGA values due to wave trapping and resonance of the sedimentary basin layers, as compared to the lower PGA values along the areas of the crystalline basement. Spectral accelerations at 10 stations show different orders of amplification over the entire period range (0.8 to 5 s) corresponding to the geological locations of these stations. The consistent decrease in amplitude with an increasing period observed at all stations indicates that it is majorly controlled by the selected source spectra. S2, S3 and S7 lie in the combined (both deep and shallow) basin area (hence, recording the maximum amplification), while S1 and S6 lie above only deepbasin area (hence, the amplification is lesser but still significant at higher periods for all three components). The rest of the stations, S4, S5, S9 and S10, are situated over the basement rocks (hence, recording the lowest value of spectral accelerations).
Our simulations focus on frequencies below 1 Hz due to high computational costs associated with sampling higher frequencies in simulations. However, this analysis remains relevant since basins, like the Kathmandu basin, often exhibit resonance at similar frequencies (Oral et al., 2022; Asimaki et al., 2017). Additionally, when dealing with higher frequencies, it becomes necessary to account for other nonlinear site effects that play a significant role in intensity variations (Semblat et al., 2005), which are not included in this analysis. More discussion on basin resonance is provided in the Supplement (Sect. S1).
Given the geometry of the basin stretched approximately north–south (NS) whilst being much more confined along east–west (EW), the amplification of both horizontal components should be theoretically contrasting. However, the periods resolved in the simulations show the intercomponent variability is still lower than the interstation variability across different geological domains (Fig. 4c). This suggests that the geometric mean of the horizontal components of PGA at each station seems to be a usable guide to explore the amplification further discussed in this study.
The pattern of higher amplification along the river basin and lower amplification along the basement area is common for PGA maps of all the earthquake scenarios (Fig. 4a). Hence while the absolute PGA is strongly dependent on the source magnitude and distance, the relative amplitude within any map is qualitatively independent of the earthquake source orientation and even magnitude. The structural similarity of PGA maps in Fig. 4a seems to indicate the potential utility of the Δ−A decomposition.
To extract this pervasive feature of relative amplification from all earthquake scenarios, we normalize and stack the PGA maps for each event. First, all PGA maps are normalized using the mean smoothearth expectation value Δ_{r}, calculated from Eq. (10). This normalization is the practical implementation of the theoretical description given in the Eq. (6), where the normalization factor is taken as the mean intensity decay in Eq. (4). Let U_{ij} be the simulated PGA at a particular site j due to earthquake i at distance r, then the normalized PGA $\left\left(\right),\widehat{{U}_{ij}}\right)$ would be
After normalization, the average PGA of the normalized maps is calculated for the number of earthquake scenarios N_{e}, as described in Eq. (7). This final, averaged PGA map is a characteristic spatial kernel for the chosen city domain and theoretically contains the average local amplification (A_{j}) at any site j for any possible earthquake regardless of source (see Fig. 5a). Here, A_{j} has the form
The calculation of A_{j} results in a mean amplification field consistent with the spatial variations observed in the simulations (Fig. 5a). Each pixel represents the mean amplification experienced at that location over all magnitudes, azimuths and directivity.
There is, of course, a dispersion of $\mathrm{ln}\left\left(\right),\widehat{{U}_{ij}}\right)$ values around this mean which is itself a spatially variable field over the domain, calculated by ${\mathit{\sigma}}_{\mathrm{ln}\left\left(\widehat{{U}_{ij}}\right\right)}$ (Fig. 5b) as
where ${\mathit{\sigma}}_{\mathrm{ln}\left(\right)open="">\left(\right)close="">\widehat{{U}_{ij}}}$ gives the variability due to various source scenarios used in the analysis and the corresponding path effects. The maximum value of ${\mathit{\sigma}}_{\mathrm{ln}\left\widehat{{U}_{ij}}\right}$ is 0.56, which is 23.8 % of the entire lnA_{j} range of 2.35 in Tomorrowville. The difference of 2.35 in maximum (lnA_{j,max}) and minimum (lnA_{j,min}) values would mean that the ratio ${A}_{j,\mathrm{max}}/{A}_{j,\mathrm{min}}$ is e^{2.35}∼ 10.48, implying an orderofmagnitude variation within Tomorrowville. Notably, the ranges of the amplification and standard deviations are of a realistic order often found in some of the extensively studied realworld settings as well, for example as shown by Day et al. (2008) in southern California.
Another approach to understanding the variability in the amplification field involves varying the number of events used to calculate lnA_{j} and examining its variability at a specific location using the events selected through a bootstrapping approach. We chose two stations from Fig. 4b, one representing an area of high amplification over the river basin, named S2, and one in low amplification over the outcropping basement, named S9 (see Fig. 5a). The number of events N_{c}, used to estimate A_{j}, is plotted against lnA_{j}, where the colour intensity represents the distribution of the iterations across the entire lnA_{j} range (Fig. 5c). For each N_{c} value, 100 random combinations of events with repetition are used for lnA_{j} calculation. The red dashes correspond to the ±1σ_{s2} and ±1σ_{s9} variability around the mean lnA_{j} value for the respective N_{c} value. The convergence of the lnA_{j} values can be observed even with as low as about seven events, with a stable ±σ_{s2} and ±σ_{s9} around the lnA_{j} values of 0.12 each. This distribution of lnA_{j} is nonoverlapping for both sites S2 and S9, which suggests that the local crustal features at both of these sites are the dominant contributor in the amplification.
The theoretical treatment described in Sect. 2 above suggests that the ground motion at a point can be decomposed into the effect of the meanfield attenuation over the wave path integrated over the crustal volume and the effect of the local velocity structure. This implies that the reversal of this process should reproduce the original PGA field. Thus if we have robust estimates of Δ and A, then we should be able to reproduce the intensity at any point using Eq. (9).
We demonstrate this process for a single earthquake, EQ13 located 30.4 km to the NW of Tomorrowville, we will show that the choice of the earthquake is not important. The simulated PGA at every point will be referred to as the true value PGA_{true} (see Fig. 6a, e). To estimate the PGA value explained in Eq. (9) for this event, referred herein as PGA_{ΔA}, we first calibrate Δ (Fig. 6b) and A (Fig. 6c) using the rest of the 39 simulated events. Δ and A are multiplied as shown in Eq. (9) to obtain PGA_{ΔA} values for this earthquake (see Fig. 6d). The difference between PGA_{ΔA} and PGA_{true} is calculated and plotted as a residual map (see Fig. 6f). The basin area shows higher negative residuals, suggesting underestimation of PGA_{ΔA} where PGA_{true} values are higher, while surrounding the basement exhibits positive values, suggesting overestimation. A graph of PGA_{ΔA} as a function of PGA_{true} is shown in Fig. 6g along with the histograms of all the grid points across Tomorrowville. There is a systematic overestimation of PGA_{ΔA} values for this particular event at the lowerPGA range, and a minor underestimation can be seen on the higherPGA side. This pattern can be attributed to the characteristic that the lnA_{j} values, which are used to calculate PGA_{ΔA}, have mean amplification values spanning a wider range compared to this specific event. The Pearson correlation coefficient (γ) between logarithms of PGA_{ΔA} and PGA_{true} is 0.98, suggesting a strong correlation between the two. The histograms presented in parallel to the axes also indicate that the distribution nature of PGA remains preserved across Tomorrowville, exhibiting a trimodal pattern in both PGA_{true} and PGA_{ΔA} (Fig. 6g). This trimodal pattern is a distinctive influence of three geological domains in the city, the deepbasin area (to the left of shallowbasin boundary), the area comprising both deep and shallow basins, and the basement region.
Finally, for each event in the suite of 40 earthquakes, the remaining 39 simulations are used to calculate Δ and A, which are multiplied to obtain PGA_{ΔA}. The results are compared with the corresponding PGA_{true} of each earthquake using the γ value and bestfitting regression line (Fig. 7a). The lowest γ value is 0.89, which suggests that the correlation is strong for all the earthquakes. In conclusion, there is a clear potential of predictability in PGA_{ΔA}, with some variability translated from different sourcespecific variability due to the heterogeneous moment distribution along the fault surface, as well as pathrelated variability due to the azimuth of sources with respect to Tomorrowville. This variability in PGA_{ΔA} is captured earlier using the ${\mathit{\sigma}}_{\mathrm{ln}\left\widehat{{U}_{ij}}\right}$ values calculated in Fig. 5b.
The impact of source orientation on the obtained γ value is illustrated by examining three parameters: the epicentral distance, back azimuth of the earthquake (bearing of the line joining hypocentre to the centre of Tomorrowville) and angle of approach (the azimuthal difference between the line connecting the hypocentre to the major fault asperity and the line connecting the hypocentre to the centre of Tomorrowville) (Fig. 7b). The back azimuth and angle of approach provide insights into the influence of the horizontally anisotropic crustal domain and directivity effects resulting from variations in fault orientation relative to Tomorrowville, respectively. γ is observed to have a positive trend with epicentral distance indicating that the earthquakes closer to Tomorrowville are poorly constrained by PGA_{ΔA} compared to the ones farther away. It can also be seen that the chosen earthquake distribution samples a wide range of values for the back azimuth and angle of approach, indicating a comprehensive representation of these factors. γ does not show any notable trend with the these two factors; hence, their impact on estimating the distribution of PGA values across Tomorrowville is not substantial.
Estimates from the United Nations Office for Disaster Risk Reduction (UNDRR) suggest that the number of people at risk from a major earthquake will increase from some 370 million in 2020 to more than 850 million by 2050 (United Nations, 2022). Due to historically unprecedented rapid urbanization, these people will be increasingly concentrated in urban centres; the same source estimates that by 2050 global urban population will increase from the current 56 % to around 68 %, with 95 % of this growth happening in the Global South. Without a concerted effort to provide decision support for risksensitive construction with a high cost–benefit ratio, ongoing urbanization in areas of high seismic hazard will increase disaster risk for millions.
That the intensity of seismic shaking varies at high spatial frequencies is graphically demonstrated by large differences in seismic damage over very short distances in areas with a uniform building code (Bielak et al., 1999; see also Asimaki et al., 2012; Dolce et al., 2003; Ohsumi et al., 2016; Sextos et al., 2018). What is less well known is the extent to which this variability is the result of differences in the earthquake source or in contrast to the rheological properties of the near surface that might impose a stable and estimable LF amplification, to first order, independent of that source. The former prioritizes forecasting likely earthquake sources in seismic hazard assessment, while the latter suggests that measuring the properties of the near surface might produce a pathway to understanding spatial patterns of seismic shaking regardless of the source. This would in turn open a path to the development of physicsbased, highresolution buildingcode classification and support evidencebased seismic urban planning policy.
Current methods for seismic hazard assessment require seismic catalogues built from the longterm deployment of large numbers of seismometers to calibrate ground motion models (Douglas and Aochi, 2008; Douglas, 2017; Douglas and Edwards, 2016). The observed variability around these models is assumed to be stochastic, and statistical methods are used to provide the moments of the emerging distributions leading to lowspatialresolution estimates of seismic hazard. Over most of the Global South such longterm data have not been collected nor is there any current appetite for deploying dense networks of seismometers required for this assessment at the resolution which would be required to guide seismicriskinformed urban planning at actionable scales.
In this study we have harnessed the potential of highresolution PB earthquake simulations to explore the extent to which seismic intensity variability might be described by nearsurface geology and explained that relative seismic intensity is independent of the earthquake source. Do some areas shake more than others, regardless of the earthquake? We exploit the certainty of a virtual world, Tomorrowville, in which the rheology, described by the geometry of the seismic velocity, is known everywhere; in which seismic sources are precisely described by kinematic models (Graves and Pitarka, 2010; Schmedes et al., 2013); and in which wave propagation is perfectly described by the wave propagation solver (SPEED) we use (Mazzieri et al., 2013). The choice of software should not lead to any notable deviation from the results obtained in this study.
The study develops a Δ−A decomposition that splits the seismic process into a meanfield attenuation model, describing the amplitude decay with source–receiver distance, and an amplification field, describing the integrated amplification of the entire wave path as experienced at each point on the surface. We have shown methods for the estimation of the Δ model and for the A field for Tomorrowville and demonstrated that their description can be used to estimate the true PGA field.
This study utilizes PB simulations in a virtual environment that shows a significant fraction of the observed variability can be explained without categorizing them as stochastic. In the real world, beyond these deterministic variations, stochastic elements of the process must be considered separately. Moreover, it becomes important to classify uncertainties as aleatory or epistemic when the real data guide the model fitting and resulting deviations (Der Kiureghian and Ditlevsen, 2009). However, in this study, PB simulation results are assumed to be devoid of any modelling uncertainties (or aleatory variability), and they are treated as reproducible true solutions in the analysis. Consequently, the deviations obtained in the results of Fig. 7a are fundamentally epistemological. The difference between the amplification map for any event and the A field that determines the value of the local PGA is precisely quantified and accessible. Investigations show that the maximum standard deviation of the A field is about 23.8 % of lnA_{j} measured across the entire area that includes the source and pathdependent variability. More importantly, analysis of the variability in the amplification value at any point indicated stable convergence from as few as seven event simulations. Furthermore, comparisons of amplifications at locations over the river basin with locations on the basement in Tomorrowville produced stable, orderofmagnitude differences in amplification which converged rapidly and which gave stable nonoverlapping amplification estimates. Of course, both the stability and the contrast in amplification are functions of the choice of velocity distribution, but the choice of the model here was developed to reflect not uncommon velocity geometry, not to accentuate amplification contrasts. We expect that the general conclusions of this work are independent of the details of the Tomorrowville velocity model.
We have not attempted to explore the variability in the amplification with the source parameters, and the initial results suggest that the influence is not likely to be strong. The main candidates expected to be dominant in the strongly anisotropic velocity model used here, source directivity and epicentral azimuth, do not make an appreciable systematic contribution to the A field. Descriptions of active fault geometry and seismotectonics of Tomorrowville could impose a source fabric introducing some systematic influence on the amplification field. Incorporation of any such influence could only constrain the variability, so the results described here might be considered a lower bound on the stability of the A field. The primary factor influencing ground motion amplification in this study is the basin geometry or buried topography, although the impact of surface topography is also anticipated to significantly affect the amplification pattern (Poursartip et al., 2020; Lee et al., 2009; GarcíaPérez et al., 2021; Geli et al., 1988). The surface topography, often rich in highresolution data, is the most straightforward to control, and it is expected to contribute to the observed variability. Future research will concentrate on investigating the influence of surface topographic features, in addition to buried topography, on the amplification phenomenon.
The reconstruction of the simulated PGA fields provided further evidence of the efficacy of the method. Using estimates of the Δ and A components from a set of 39 simulations (out of 40) provided strong correlations between true and inverted PGA fields for the 40th. Further, in keeping with the observation of nonoverlapping amplification values for basement and basin locations, places with highintensity shaking were broadly consistently high intensity for all events, whereas locations experiencing lowintensity shaking were also consistent across all events.
The results are suggestive of an underlying physical process in which smallscale LF relative shaking intensity is controlled more by local geology than by source process. Given the description of the relevant fields through simulations, each taking approximately a day on a commonly available computer cluster (see Table S3 for simulation parameters and runtime estimates), it is feasible to estimate the entire PGA field (PGA_{ΔA}) for an event of a specific magnitude and location in milliseconds of computing time. At a minimum, this provides a workflow through which normal probabilistic seismic hazard assessments that require estimates of PGA for thousands of events at each location can benefit from the advances in PB simulations without the massive compute overhead that makes these computations unfeasible at present.
The stability of the relative amplification field together with the stable, orderofmagnitude difference in PGA across the surface of Tomorrowville demonstrated in this study points to methods for highresolution seismic hazard estimation based on understanding the static properties of the near surface, rather than on the unpredictable properties of future earthquakes. The challenge becomes a problem of measurement, rather than forecasting. There remains the critical problem either of the elucidation of the velocity structure of the near surface (Sebastiano et al., 2019), so the Δ and A fields might be estimated through simulation as in this paper, or of the direct estimation of the field by measurement of the intensity of shaking at high resolution in the area of interest. To clarify again, this study explores only LF nearsurface effects arising from the presence of complex sedimentary basins and shows their contribution in shortscale variability in amplification. It is noteworthy that these LF effects are additional to the site effects related to very nearsurface (decametre) depths, which include nonlinear soil responses and other highspatialfrequency velocity variations, all of which can lead to intricate outcomes (Taborda et al., 2012). Consequently, for applications like enhancing microzonation maps, it is imperative to merge this analysis with elements accounting for HF variability.
In conclusion, rapid urban expansion in areas of poor historical instrumentation leaves significant gaps in data for seismic hazard assessment. Furthermore, current methods both require decadelong deployment of dense seismic networks in the area of nearfuture urban development and fail to provide highresolution assessments that identify areas of strong and weak shaking that could underpin seismiccode classification with a high cost–benefit ratio. The potential of PB simulations has prompted the evaluation of the seismicwave field across areas of nearfuture development. The results suggest methods to allow for the rapid, highresolution assessment of geological structure that could lead to risk assessment at unprecedented resolution.
The data used in this research are mainly the simulation outputs, which are extensive in scale. The critical information regarding the crustal domain; earthquake hypocentre; and PGA data, which are pivotal for generating the majority of the paper's results, can be found in the Supplement. For more detailed information on earthquake moment distribution, we encourage readers to refer to Jenkins et al. (2023). The software used to run the simulation is an opensource package, SPEED (https://speed.mox.polimi.it/download/, Mazzieri et al., 2013). The data analysis and processing are done using Python, and the code is available at https://github.com/himansh78/GroundMotionCalc.git (last access: 10 October 2024; https://doi.org/10.5281/zenodo.13884828, Agrawal, 2024).
The supplement related to this article is available online at: https://doi.org/10.5194/nhess2435192024supplement.
HA and JM contributed to the study conception and design. Material preparation and data analysis were performed by HA. The first draft of the paper was prepared by HA, including all the figures and text. The text was further reviewed and improved with the supervision of JM. JM passed away before submission.
The contact author has declared that none of the authors has any competing interests.
Publisher’s note: Copernicus Publications remains neutral with regard to jurisdictional claims made in the text, published maps, institutional affiliations, or any other geographical representation in this paper. While Copernicus Publications makes every effort to include appropriate place names, the final responsibility lies with the authors.
John McCloskey is listed as a coauthor in recognition of his significant contributions. Unfortunately, he passed away after the paper was ready for submission, and we deeply mourn his loss.
The authors acknowledge initial discussions and simulations obtained with the prompt support and guidance from Karim Tarbali, a former postdoctoral research associate (PDRA) at the University of Edinburgh. We thank Gemma Cremen, Chris J. Bean, Mark Naylor, Ian Main, Karen Lythgoe and two anonymous reviewers for providing constructive feedback and guidance that improved the paper.
This research is a part of the wider PhD project “Physicsbased ground motion simulations for seismic risk assessment and uncertainty quantification in rapidly urbanising environments”. The PhD student is funded by the University of Edinburgh School of GeoSciences. This research project is also supported by the Tomorrow's Cities hub (UKRI GCRF fund under grant no. NE/S009000/1).
This paper was edited by Solmaz Mohadjer and reviewed by two anonymous referees.
Abrahamson, N. A., Silva, W. J., and Kamai, R.: Summary of the ASK14 ground motion relation for active crustal regions, Earthq. Spectra, 30, 1025–1055, https://doi.org/10.1193/070913EQS198M, 2014.
Agrawal, H.: Supplementary code to “Estimating ground motion intensities using simulationbased estimates of local crustal seismic response” (v1.0.0), Zenodo [code], https://doi.org/10.5281/zenodo.13884828, 2024.
Aki, K.: Local site effects on weak and strong ground motion, Tectonophysics, 218, 93–111, https://doi.org/10.1016/00401951(93)90262I, 1993.
Ancheta, T. D., Darragh, R. B., Stewart, J. P., Seyhan, E., Silva, W. J., Chiou, B. S. J., Wooddell, K. E., Graves, R. W., Kottke, A. R., Boore, D. M., Kishida, T., and Donahue, J. L.: NGAWest2 database, Earthq. Spectra, 30, 989–1005, https://doi.org/10.1193/070913EQS197M, 2014.
Anderson, J. G. and Brune, J. N.: Probabilistic seismic hazard analysis without the ergodic assumption, Seismol. Res. Lett., 70, 19–28, https://doi.org/10.1785/gssrl.70.1.19, 1999.
Ansal, A., Kurtuluş, A., and Tönük, G.: Seismic microzonation and earthquake damage scenarios for urban areas, Soil Dyn. Earthq. Eng., 30, 1319–1328, https://doi.org/10.1016/j.soildyn.2010.06.004, 2010.
Asimaki, D., Ledezma, C., Montalva, G. A., Tassara, A., Mylonakis, G., and Boroschek, R.: Site effects and damage patterns, Earthq. Spectra, 28, S55–S74, https://doi.org/10.1193/1.4000029, 2012.
Asimaki, D., Mohammadi, K., Mason, H. B., Adams, R. K., Rajaure, S., and Khadka, D.: Observations and Simulations of Basin Effects in the Kathmandu Valley during the 2015 Gorkha, Nepal, Earthquake Sequence, Earthq. Spectra, 33, S35–S53, https://doi.org/10.1193/013117EQS022M, 2017.
Atkinson, G. M. and Boore, D. M.: Earthquake GroundMotion Prediction Equations for Eastern North America, B. Seismol. Soc. Am., 96, 2181–2205, https://doi.org/10.1785/0120050245, 2006.
Baker, J. W., Bradley, B. A., and Stafford, P. J.: Probabilistic seismic hazard and risk analysis, Cambridge University Press, Cambridge, https://doi.org/10.1017/9781108425056, 2021.
Bazzurro, P. and Cornell, C. A.: Nonlinear soilsite effects in probabilistic seismichazard analysis, B. Seismol. Soc. Am., 94, 2110–2123, https://doi.org/10.1785/0120030216, 2004.
Bielak, J., Xu, J., and Ghattas, O.: Earthquake Ground Motion and Structural Response in Alluvial Valleys, J. Geotech. Geoenviron., 125, 413–423, https://doi.org/10.1061/(ASCE)10900241(1999)125:5(413), 1999.
Borcherdt, R. D. and Glassmoyer, G.: On the characteristics of local geology and their influence on ground motions generated by the Loma Prieta earthquake in the San Franciso Bay region, California, B. Seismol. Soc. Am., 82, 603–641, https://doi.org/10.1785/bssa0820020603, 1992.
Bradley, B. A.: A framework for validation of seismic response analyses using seismometer array recordings, Soil Dyn. Earthq. Eng., 31, 512–520, https://doi.org/10.1016/j.soildyn.2010.11.008, 2011.
Bradley, B. A.: Systematic ground motion observations in the Canterbury earthquakes and regionspecific nonergodic empirical ground motion modeling, Earthq. Spectra, 31, 1735–1761, https://doi.org/10.1193/053013EQS137M, 2015.
Bradley, B. A.: Ongoing challenges in physicsbased ground motion prediction and insights from the 2010–2011 Canterbury and 2016 Kaikoura, New Zealand earthquakes, Soil Dyn. Earthq. Eng., 124, 354–364, https://doi.org/10.1016/j.soildyn.2018.04.042, 2019.
Brocher, T. M.: Empirical relations between elastic wavespeeds and density in the Earth's crust, B. Seismol. Soc. Am., 95, 2081–2092, https://doi.org/10.1785/0120050077, 2005.
Brocher, T. M.: Compressional and shearwave velocity versus depth relations for common rock types in northern California, B. Seismol. Soc. Am., 98, 950–968, https://doi.org/10.1785/0120060403, 2008.
Cadet, H., Macau, A., Benjumea, B., Bellmunt, F., and Figueras, S.: From ambient noise recordings to site effect assessment: The case study of Barcelona microzonation, Soil Dyn. Earthq. Eng., 31, 271–281, https://doi.org/10.1016/j.soildyn.2010.07.005, 2011.
Campbell, K. W. and Bozorgnia, Y.: NGAWest2 ground motion model for the average horizontal components of PGA, PGV, and 5 % damped linear acceleration response spectra, Earthq. Spectra, 30, 1087–1114, https://doi.org/10.1193/062913EQS175M, 2014.
Castellaro, S. and Musinu, G.: Resonance versus Shape of Sedimentary Basins, B. Seismol. Soc. Am., 113, 745–761, https://doi.org/10.1785/0120210277, 2023.
Castellaro, S., Mulargia, F., and Rossi, P. L.: Vs30: Proxy for seismic amplification?, Seismol. Res. Lett., 79, 540–543, https://doi.org/10.1785/gssrl.79.4.540, 2008.
Chiou, B. S. J. and Youngs, R. R.: Update of the Chiou and Youngs NGA model for the average horizontal component of peak ground motion and response spectra, Earthq. Spectra, 30, 1117–1153, https://doi.org/10.1193/072813EQS219M, 2014.
Cramer, C. H.: Sitespecific seismichazard analysis that is completely probabilistic, B. Seismol. Soc. Am., 93, 1841–1846, https://doi.org/10.1785/0120020206, 2003.
Cremen, G., Galasso, C., McCloskey, J., Barcena, A., Creed, M., Filippi, M. E., Gentile, R., Jenkins, L. T., Kalaycioglu, M., Mentese, E. Y., Muthusamy, M., Tarbali, K., and Trogrlić, R. Š.: A stateoftheart decisionsupport environment for risksensitive and propoor urban planning and design in Tomorrow's cities, Int. J. Disast. Risk Re., 85, 103400, https://doi.org/10.1016/j.ijdrr.2022.103400, 2023.
Day, S. M., Graves, R., Bielak, J., Dreger, D., Larsen, S., Olsen, K. B., Pitarka, A., and RamirezGuzman, L.: Model for Basin Effects on LongPeriod Response Spectra in Southern California:, Earthq. Spectra, 24, 257–277, https://doi.org/10.1193/1.2857545, 2008.
De Hoop, A. T.: Representation theorems for the displacement in an elastic solid and their application to elastodynamic diffraction theory, Technische Hogeschoo, Delft, Netherlands, https://repository.tudelft.nl/record/uuid:4678f6dc23c24705b6b8d22bbc4921e1 (last access: 3 October 2024), 1958.
Der Kiureghian, A. and Ditlevsen, O.: Aleatory or epistemic? Does it matter?, Struct. Saf., 31, 105–112, https://doi.org/10.1016/j.strusafe.2008.06.020, 2009.
Dolce, M., Masi, A., Marino, M., and Vona, M.: Earthquake damage scenarios of the building stock of Potenza (Southern Italy) including site effects, B. Earthq. Eng., 1, 115–140, https://doi.org/10.1023/A:1024809511362, 2003.
Douglas, J.: Ground motion prediction equations 1964–2019 (December 2019), SED Rep. SED/ENSI/R/01/20140911, 651 pp., https://www.linkedin.com/posts/johndouglasglasgow_gmpegroundmotionmodelseismichazardactivity7233422061291929600Rp5m?utm_source=share&utm_medium=member_desktop (last access: 3 October 2024), 2017.
Douglas, J. and Aochi, H.: A survey of techniques for predicting earthquake ground motions for engineering purposes, Surv. Geophys., 29, 187–220, https://doi.org/10.1007/s107120089046y, 2008.
Douglas, J. and Edwards, B.: Recent and future developments in earthquake ground motion estimation, EarthSci. Rev., 160, 203–219, https://doi.org/10.1016/j.earscirev.2016.07.005, 2016.
Foti, S., Aimar, M., Ciancimino, A., and Passeri, F.: Recent developments in seismic site response evaluation and microzonation, in: Proceedings of the 17th European Conference on Soil Mechanics and Geotechnical Engineering, Reykjavik, Iceland, 1–6 September 2019, ECSMGE 2019, https://doi.org/10.32075/17ECSMGE20191117, 2019.
Frankel, A.: Threedimensional simulations of ground motions in the San Bernardino Valley, California, for hypothetical earthquakes on the San Andreas Fault, B. Seismol. Soc. Am., 83, 1020–1041, https://doi.org/10.1785/BSSA0830041020, 1993.
Freddi, F., Galasso, C., Cremen, G., Dall'Asta, A., Di Sarno, L., Giaralis, A., GutiérrezUrzúa, F., MálagaChuquitaype, C., Mitoulis, S. A., Petrone, C., Sextos, A., Sousa, L., Tarbali, K., Tubaldi, E., Wardman, J., and Woo, G.: Innovations in earthquake risk reduction for resilience: Recent advances and challenges, Int. J. Disast. Risk Re., 60, 102267, https://doi.org/10.1016/j.ijdrr.2021.102267, 2021.
GarcíaPérez, T., Ferreira, A. M. G., Yáñez, G., Iturrieta, P., and Cembrano, J.: Effects of topography and basins on seismic wave amplification: The Northern Chile coastal cliff and intramountainous basins, Geophys. J. Int., 227, 1143–1167, https://doi.org/10.1093/gji/ggab259, 2021.
Geli, L., Bard, P. Y., and Jullien, B.: The effect of topography on earthquake ground motion: A review and new results, B. Seismol. Soc. Am., 78, 42–63, https://doi.org/10.1785/BSSA0780010042, 1988.
Gentile, R., Cremen, G., Galasso, C., Jenkins, L. T., Manandhar, V., Mentese, E. Y., Guragain, R., and McCloskey, J.: Scoring, selecting, and developing physical impact models for multi hazard risk assessment, Int. J. Disast. Risk Re., 82, 103365, https://doi.org/10.1016/j.ijdrr.2022.103365, 2022.
Graves, R. W. and Pitarka, A.: Broadband groundmotion simulation using a hybrid approach, B. Seismol. Soc. Am., 100, 2095–2123, https://doi.org/10.1785/0120100057, 2010.
Graves, R. W., Pitarka, A., and Somerville, P. G.: Groundmotion amplification in the Santa Monica area: Effects of shallow basinedge structure, B. Seismol. Soc. Am., 88, 1224–1242, https://doi.org/10.1785/bssa0880051224, 1998.
Hough, S. E. and Anderson, J. G.: Highfrequency Spectra Observed at Anza, California: Implications for Q Structure, B. Seismol. Soc. Am., 78, 692–707, https://doi.org/10.1785/BSSA0780020692, 1988.
Hough, S. E., Martin, S. S., Gahalaut, V., Joshi, A., Landes, M., and Bossu, R.: A comparison of observed and predicted ground motions from the 2015 M_{W} 7.8 Gorkha, Nepal, earthquake, Nat. Hazards, 84, 1661–1684, https://doi.org/10.1007/s1106901625058, 2016.
Jenkins, L. T., Creed, M. J., Tarbali, K., Muthusamy, M., Trogrlić, R. Š., Phillips, J. C., Watson, C. S., Sinclair, H. D., Galasso, C., and McCloskey, J.: Physicsbased simulations of multiple natural hazards for risksensitive planning and decisionmaking in expanding urban regions, Int. J. Disast. Risk Re., 84, 103338, https://doi.org/10.1016/j.ijdrr.2022.103338, 2023.
Kaklamanos, J., Bradley, B. A., Thompson, E. M., and Baise, L. G.: Critical parameters affecting bias and variability in siteresponse analyses using KiKnet downhole array data, B. Seismol. Soc. Am., 103, 1733–1749, https://doi.org/10.1785/0120120166, 2013.
Kamai, R., Abrahamson, N. A., and Silva, W. J.: VS30 in the NGA GMPEs: Regional differences and suggested practice, Earthq. Spectra, 32, 2083–2108, https://doi.org/10.1193/072615EQS121M, 2016.
Knopoff, L.: Diffraction of Elastic Waves, J. Acoust. Soc. Am., 28, 217–229, https://doi.org/10.1121/1.1908247, 1956.
Kramer, S. L.: Geotechnical Earthquake Engineering, Pearson PrenticeHall, Upper Saddle River, NJ, USA, https://faculty.washington.edu/kramer/GEEbook.pdf (last access: 3 October 2024), 1996.
Kramer, S. L. and Mitchell, R. A.: Ground motion intensity measures for liquefaction hazard evaluation, Earthq. Spectra, 22, 413–438, https://doi.org/10.1193/1.2194970, 2006.
Kuehn, N. M., Abrahamson, N. A., and Walling, M. A.: Incorporating Nonergodic Path Effects into the NGAWest2 GroundMotion Prediction Equations, B. Seismol. Soc. Am., 109, 575–585, https://doi.org/10.1785/0120180260, 2019.
Landwehr, N., Kuehn, N. M., Scheffer, T., and Abrahamson, N.: A Nonergodic GroundMotion Model for California with Spatially Varying Coefficients, B. Seismol. Soc. Am., 6, 2574–2583, https://doi.org/10.1785/0120160118, 2016.
Lee, S. J., Komatitsch, D., Huang, B. S., and Tromp, J.: Effects of topography on seismicwave propagation: An example from Northern Taiwan, B. Seismol. Soc. Am., 99, 314–325, https://doi.org/10.1785/0120080020, 2009.
Liu, P., Archuleta, R. J., and Hartzell, S. H.: Prediction of broadband groundmotion time histories: Hybrid low/highfrequency method with correlated random source parameters, B. Seismol. Soc. Am., 96, 2118–2130, https://doi.org/10.1785/0120060036, 2006.
Marafi, N. A., Eberhard, M. O., Berman, J. W., Wirth, E. A., and Frankel, A. D.: Effects of deep basins on structural collapse during large subduction earthquakes, Earthq. Spectra, 33, 963–997, https://doi.org/10.1193/071916EQS114M, 2017.
Maufroy, E., CruzAtienza, V. M., and Gaffet, S.: A robust method for assessing 3D topographic site effects: A case study at the LSBB underground laboratory, France, Earthq. Spectra, 28, 1097–1115, https://doi.org/10.1193/1.4000050, 2012.
Mazzieri, I., Stupazzini, M., Guidotti, R., and Smerzini, C.: SPEED: SPectral elements in elastodynamics with discontinuous galerkin: A nonconforming approach for 3d multiscale problems, Int. J. Numer. Meth. Eng., 95, 991–1010, https://doi.org/10.1002/nme.4532, 2013 (code available at: https://speed.mox.polimi.it/download/, last access: 3 October 2024).
McGuire, R. K.: Probabilistic seismic hazard analysis: Early history, Earthq. Eng. Struct. D., 37, 329–338, https://doi.org/10.1002/eqe.765, 2008.
Menteşe, E. Y., Cremen, G., Gentile, R., Galasso, C., Filippi, E. M., and McCloskey, J.: Future exposure modelling for riskinformed decision making in urban planning, Int. J. Disast. Risk Re., 90, 103651, https://doi.org/10.1016/j.ijdrr.2023.103651, 2023.
Mucciarelli, M. and Gallipoli, M. R.: Comparison between V_{s30} and other estimates of site amplification in Italy, in: First European Conference on Earthquake Engineering and Seismology, Geneva, Switzerland, 3–8 September 2006, 270, http://hdl.handle.net/2122/1945 (last access: 3 October 2024), 2006.
Nath, S. K. and Thingbaijam, K. K. S.: Peak ground motion predictions in India: an appraisal for rock sites, J. Seismol., 15, 295–315, https://doi.org/10.1007/s1095001092245, 2011.
Ohsumi, T., Mukai, Y., and Fujitani, H.: Investigation of Damage in and Around Kathmandu Valley Related to the 2015 Gorkha, Nepal Earthquake and Beyond, Geotech. Geol. Eng., 34, 1223–1245, https://doi.org/10.1007/s1070601600239, 2016.
Oral, E., Ayoubi, P., Ampuero, J. P., Asimaki, D., and Bonilla, L. F.: Kathmandu Basin as a local modulator of seismic waves: 2D modelling of nonlinear site response under obliquely incident waves, Geophys. J. Int., 231, 1996–2008, https://doi.org/10.1093/gji/ggac302, 2022.
Paolucci, R., Mazzieri, I., Smerzini, C., and Stupazzini, M.: PhysicsBased Earthquake Ground Shaking Scenarios in Large Urban Areas, in: Geotechnical, Geological and Earthquake Engineering, vol. 34, Kluwer Academic Publishers, 331–359, https://doi.org/10.1007/9783319071183_10, 2014.
Pilz, M., Parolai, S., Stupazzini, M., Paolucci, R., and Zschau, J.: Modelling basin effects on earthquake ground motion in the Santiago de Chile basin by a spectral element code, Geophys. J. Int., 187, 929–945, https://doi.org/10.1111/j.1365246X.2011.05183.x, 2011.
Pitilakis, K., Riga, E., Anastasiadis, A., Fotopoulou, S., and Karafagka, S.: Towards the revision of EC8: Proposal for an alternative site classification scheme and associated intensity dependent spectral amplification factors, Soil Dyn. Earthq. Eng., 126, 105137, https://doi.org/10.1016/j.soildyn.2018.03.030, 2019.
Poursartip, B., Fathi, A., and Tassoulas, J. L.: Largescale simulation of seismic wave motion: A review, Soil Dyn. Earthq. Eng., 129, 105909, https://doi.org/10.1016/j.soildyn.2019.105909, 2020.
RodriguezMarek, A., Rathje, E. M., Bommer, J. J., Scherbaum, F., and Stafford, P. J.: Application of singlestation sigma and siteresponse characterization in a probabilistic seismichazard analysis for a new nuclear site, B. Seismol. Soc. Am., 104, 1601–1619, https://doi.org/10.1785/0120130196, 2014.
Schmedes, J., Archuleta, R. J., and Lavalĺee, D.: A kinematic rupture model generator incorporating spatial interdependency of earthquake source parameters, Geophys. J. Int., 192, 1116–1131, https://doi.org/10.1093/gji/ggs021, 2013.
Sebastiano, D., Francesco, P., Salvatore, M., Roberto, I., Antonella, P., Giuseppe, L., Pauline, G., and Daniela, F.: Ambient noise techniques to study nearsurface in particular geological conditions: A brief review, in: Innovation in NearSurface Geophysics: Instrumentation, Application, and Data Processing Methods, Elsevier Inc., 419–460, https://doi.org/10.1016/B9780128124291.00012X, 2019.
Semblat, J. F., Kham, M., Parara, E., Bard, P. Y., Pitilakis, K., Makra, K., and Raptakis, D.: Seismic wave amplification: Basin geometry vs soil layering, Soil Dyn. Earthq. Eng., 25, 529–538, https://doi.org/10.1016/j.soildyn.2004.11.003, 2005.
Sextos, A., De Risi, R., Pagliaroli, A., Foti, S., Passeri, F., Ausilio, E., Cairo, R., Capatti, M. C., Chiabrando, F., Chiaradonna, A., Dashti, S., De Silva, F., Dezi, F., Durante, M. G., Giallini, S., Lanzo, G., Sica, S., Simonelli, A. L., and Zimmaro, P.: Local site effects and incremental damage of buildings during the 2016 Central Italy Earthquake sequence, Earthq. Spectra, 34, 1639–1669, https://doi.org/10.1193/100317EQS194M, 2018.
Shi, J. and Asimaki, D.: From Stiffness to Strength: Formulation and Validation of a Hybrid Hyperbolic Nonlinear Soil Model for SiteResponse Analyses, B. Seismol. Soc. Am., 107, 1336–1355, https://doi.org/10.1785/0120150287, 2017.
Smerzini, C. and Villani, M.: Broadband numerical simulations in complex nearfield geological configurations: The case of the 2009 M_{w} 6.3 L'Aquila earthquake, B. Seismol. Soc. Am., 102, 2436–2451, https://doi.org/10.1785/0120120002, 2012.
Smerzini, C., Paolucci, R., and Stupazzini, M.: Comparison of 3D, 2D and 1D numerical approaches to predict long period earthquake ground motion in the Gubbio plain, Central Italy, B. Earthq. Eng., 9, 2007–2029, https://doi.org/10.1007/s1051801192898, 2011.
Spudich, P., Bayless, J. R., Baker, J., Chiou, B. S. J., Rowshandel, B., Shahi, S., and Somerville, P.: Final Report of the NGAWest2 Directivity Working Group, Pacific Engineering Research Center Report, California, https://researchers.mq.edu.au/en/publications/finalreportofthengawest2directivityworkinggroup (last access: 3 October 2024), 2013.
Stewart, J. P., Afshari, K., and Goulet, C. A.: Nonergodic site response in seismic hazard analysis, Earthq. Spectra, 33, 1385–1414, https://doi.org/10.1193/081716EQS135M, 2017.
Stirling, M. W.: The Continued Utility of Probabilistic SeismicHazard Assessment, in: Earthquake Hazard, Risk and Disasters, Elsevier Inc., 359–376, https://doi.org/10.1016/B9780123948489.000134, 2014.
Stirling, M. W., McVerry, G., Gerstenberger, M., Litchfield, N., Van Dissen, R., Berryman, K., Barnes, P., Wallace, L., Villamor, P., Langridge, R., Lamarche, G., Nodder, S., Reyners, M., Bradley, B., Rhoades, D., Smith, W., Nicol, A., Pettinga, J., Clark, K., and Jacobs, K.: National seismic hazard model for New Zealand: 2010 update, B. Seismol. Soc. Am., 102, 1514–1542, https://doi.org/10.1785/0120110170, 2012.
Taborda, R. and Roten, D.: PhysicsBased GroundMotion Simulation, in: Encyclopedia of Earthquake Engineering, 33 pp., https://doi.org/10.1007/9783642361975, 2014.
Taborda, R., Bielak, J., and Restrepo, D.: Earthquake GroundMotion Simulation including Nonlinear Soil Effects under Idealized Conditions with Application to Two Case Studies, Seismol. Res. Lett., 83, 1047–1060, https://doi.org/10.1785/0220120079, 2012.
Torre, C. A. de la, Bradley, B. A., and Lee, R. L.: Modeling nonlinear site effects in physicsbased ground motion simulations of the 2010–2011 Canterbury earthquake sequence, Earthq. Spectra, 36, 856–879, https://doi.org/10.1177/8755293019891729, 2020.
Tsai, C. C., Kishida, T., and Lin, W. C.: Adjustment of site factors for basin effects from site response analysis and deep downhole array measurements in Taipei, Eng. Geol., 285, 106071, https://doi.org/10.1016/j.enggeo.2021.106071, 2021.
UNISDR: Sendai framework for disaster risk reduction 20152030, Aust. J. Emerg. Manag., 30, 9–10, 2015.
United Nations: Envisaging the Future of Cities. World Cities report 2022, UN, 422 pp., https://unhabitat.org/sites/default/files/2022/06/wcr_2022.pdf (last access: 3 October 2024), 2022.
Wang, C., Cremen, G., Gentile, R., and Galasso, C.: Design and assessment of propoor financial soft policies for expanding cities, Int. J. Disast. Risk Re., 85, 103500, https://doi.org/10.1016/j.ijdrr.2022.103500, 2023.
Wang, G., Du, C., Huang, D., Jin, F., Koo, R. C. H., and Kwan, J. S. H.: Parametric models for 3D topographic amplification of ground motions considering subsurface soils, Soil Dyn. Earthq. Eng., 115, 41–54, https://doi.org/10.1016/j.soildyn.2018.07.018, 2018.
Wells, D. L. and Coppersmith, K. J.: New empirical relationships among magnitude, rupture length, rupture width, rupture area, and surface displacement, B. Seismol. Soc. Am., 84, 974–1002, https://doi.org/10.1785/bssa0840040974, 1994.
Yomogida, K. and Etgen, J. T.: 3D wave propagation in the Los Angeles Basin for the WhittierNarrows earthquake, B. Seismol. Soc. Am., 83, 1325–1344, https://doi.org/10.1086/622062, 1993.
Zhu, C., Thambiratnam, D., and Gallage, C.: Statistical analysis of the additional amplification in deep basins relative to the 1D approach, Soil Dyn. Earthq. Eng., 104, 296–306, https://doi.org/10.1016/j.soildyn.2017.09.003, 2018.
 Abstract
 Key points
 Introduction
 Theoretical considerations
 Defining domain and source scenarios for simulations
 Estimation of Δ and A for Tomorrowville
 Estimation of PGA using Δ and A for 40 earthquakes
 Discussion and summary
 Code and data availability
 Author contributions
 Competing interests
 Disclaimer
 Acknowledgements
 Financial support
 Review statement
 References
 Supplement
 Abstract
 Key points
 Introduction
 Theoretical considerations
 Defining domain and source scenarios for simulations
 Estimation of Δ and A for Tomorrowville
 Estimation of PGA using Δ and A for 40 earthquakes
 Discussion and summary
 Code and data availability
 Author contributions
 Competing interests
 Disclaimer
 Acknowledgements
 Financial support
 Review statement
 References
 Supplement