Investigation of topographic site effects using 3D waveform modelling: amplification, polarization and torsional motions in the case study of Arquata del Tronto (Italy)

The combined effect of topography and near-surface heterogeneities on the seismic response is hardly predictable and may lead to an aggravation of the ground motion. We apply physics-based numerical simulations of 3D seismic wave propagation to highlight these effects in the case study of Arquata del Tronto, a municipality in the Apennines that includes a historical village on a hill and a hamlet on the flat terrain of an alluvial basin. The two hamlets suffered different damage during the 2016 seismic sequence in Central Italy. We analyze the linear visco-elastic seismic response for vertically incident plane waves in terms of spectral amplification, polarization and induced torsional motion within the frequency band 1–8 Hz over a 1 km2 square area, with spatial resolution 25 m. To discern the effects of topography from those of the sub-surface structure we iterate the numerical simulations for three different versions of the sub-surface model: one homogeneous, one with a surficial weathering layer and a soil basin and one with a complex internal setting. The numerical results confirm the correlation between topographic curvature and amplification and support a correlation between the induced torsional motion and the topographic slope. On the other hand we find that polarization does not necessarily imply ground motion amplification. In the frequency band above 4 Hz the topography-related effects are mainly aggravated by the presence of the weathering layer, even though they do not exceed the soil-related effects in the flat-topography basin. The geological setting below the weathering layer plays a recognizable role in the topography-related site response only for frequencies below 4 Hz.


Introduction
The estimation of ground motion amplification due to site effects plays a fundamental role in the efforts devoted to seismic-risk mitigation (Amanti et al. 2020). Thanks to ground motion recordings and earthquake damage reports from sites characterized by topography, it has long been known that topographic features are an indicator of site effects, but theoretical explanations based on the relief's geometrical characteristics usually underestimated the observed amplification values (Geli et al. 1988). The discrepancies between predictions and observations can be ascribed to the absence of adequate reference site in many experimental estimates (Pedersen et al. 1994;Wood and Cox 2016) on one hand, and on the other hand, to the adoption of oversimplified models in theoretical estimates. The oversimplification that may cause underestimation of theoretically predicted amplifications typically consists of inadequate two-dimensional (2D) approximation of prominently three-dimensional (3D) topographic features (Paolucci 2002;Luo et al. 2020) or in neglecting impedance contrasts in the subsurface structure that possibly concur with topography in forming amplification effects (Graizer 2009;Assimaki and Jeong 2013;Hailemikael et al. 2016, to cite a few).
There is today a general agreement that topography implies a number of effects, discussed in the overview provided by Massa et al. (2014) and briefly summarized in the points below: • Frequency dependent amplification occurs at the top of the topographic high in the wavelength range comparable to the width of the morphological relief (peak frequency decreases with increasing relief width). • Amplification is proportional to the shape ratio (height vs. base half-width) of the relief. • Amplification is larger for horizontal components, in particular for the component transverse to the direction of elongation of the relief. • Topographic amplification is dependent on the source-to-site direction.
• Topographic effects combined with the effects of subsoil heterogeneities, give rise to atypical topographic effects.
In order to account for possible topography-related amplification effects, a number of seismic codes, as Eurocode 8 (CEN 2004) and NTC18 in Italy (Ministero delle Infrastrutture e dei Trasporti 2018), have introduced topographic amplification coefficients for simplified morphological conditions. Recently the efforts are focused on the development of more sophisticated estimators, based on the availability of digital elevation models (DEM) and equivalent uniform rock stiffness estimates (Grelle et al. 2018). Zhou et al. (2020) adopt the back-propagation neural network technique for the derivation of a ground motion amplification model based on topographic geometrical features estimated from the DEM. Since they are based on surface topography alone, these approaches neglect possible interactions between topography and subsurface complexities on wave propagation. According to Burjánek et al. (2014), the observed amplification at sites with comparable topography but with different geology may differ for a factor which is significantly larger than the expected ground motion variability and studies based exclusively on the terrain topography have almost no chance to capture the site effects correctly.
In the present work we address the problem of possible ground motion aggravation due to the complex interaction between a pronounced topographic relief and a heterogeneous distribution of underground seismic velocities. Given that the distinction of topographic effects from the stratigraphic ones in the seismic records results impractical, we adopted an approach that fully relies on numerical simulations to separate the two (Ashford et al. 1997). Three-dimensional (3D) physics-based numerical simulations of seismic wave propagation can provide a detailed and accurate characterization of the seismic ground motion and allow us to explore its spatial heterogeneity. This kind of numerical simulations have proven to be essential for the explanation of anomalies in the observed ground motion during past earthquakes (e. g. Paolucci et al. 2015;Cruz-Atienza et al. 2016;Klin et al. 2019), also in cases in which topographic effects are involved (e.g. Puglia et al. 2013;Luo et al. 2020). Moreover, numerical simulations allow the quantification of potential amplification effects at sites where ground motion data is lacking (e. g. Fehr et al. 2019).
In the present work we analyze the site response not only in the conventional terms of spectral amplification but also in terms of induced polarization and induced rotations in the horizontal components of ground motion.
According to several works (e. g. Burjánek et al. 2014;Massa et al. 2014;Pischiutta et al. 2018), the polarization of the site response can be viewed as an indicator of topography-induced strong site effects, even though the direction of polarization is not strictly related to the main directions of the topographic features. Given that the polarization of the site response can be identified by low-cost single-station analyses, deepening this topic could provide useful indications for seismic microzonation studies.
On the other hand, ground motion rotations are rarely considered in site effects studies, due to difficulties involved in measuring rotational motions and strains and because of a widespread preconception in the seismological community that rotational motions are insignificant (Lee et al. 2009). In the earthquake engineering community instead, the ground rotational motion gained interest (Lee and Trifunac 1985;Castellani et al. 2012), since the travelling wave effect was recognized by Newmark (1969) and it is today accounted for in building codes under the term accidental eccentricity (Anagnostopoulos et al. 2015). While the research regarding the relevance of ground motion rotations is still ongoing (Guidotti et al. 2018), in the present work we investigate the possibility of topography-induced rotations around the vertical axis (torsional motion). Topography-related torsional motions were observed-without a quantitative analysis-by Stupazzini et al. (2009) in their 3D numerical simulations of earthquake ground motion in the Grenoble valley and further analyses were addressed in Smerzini et al. (2009). In the present work we set up a simple predictive model which relates torsional motions and the topographic slope in homogeneous media and we evaluate the possible variations due to near surface heterogeneities.
We elected as case study for the aforementioned investigations the 1 km 2 square area centered in the historical center of Arquata del Tronto, a village in the Central Apennines located on an elongated WNW-ESE ridge at an elevation of about 100 m above the underlying valley. The village suffered extensive damage and collapses due to the 2016 August 24th Mw 6.0 Amatrice earthquake and the later shocks and the non-uniform damage distribution indicated possible topography-related site effects (Galli et al. 2016;Lanzo et al. 2019). An evaluation of site effects in Arquata del Tronto and Borgo from seismological data has been performed by Laurenzano et al. (2019), who exploited earthquake recordings acquired on temporary seismic stations deployed after the first mainshock of the 2016 seismic sequence. Their study evidenced a substantial difference between the frequency-dependent site amplification function at a location on the top of the ridge and in the valley below. A directional effect that could be associated with the ridge orientation has been pointed out, as well. In order to quantitatively explain the observed site response characteristics, Giallini et al. (2020) applied 2D numerical simulations using subsoil models deduced from the available geognostic, geophysical and geomechanical surveys and evidenced the coupling effect of stratigraphy and topography associated with the presence of a weathered portion of the rock mass and the alternation of highly dipping rocky materials. They admitted that 2D analyses cannot fully capture the behavior of such an asymmetric ridge and suggested tridimensional numerical analyses for this site. This origin of the site effects was finally verified by Primofiore et al. (2020), who set up an updated 3D digital model of the area and evaluated the resulting ground motion by means of 3D numerical simulations of seismic wave propagation.
The present work can be therefore considered as a follow up of the analyses by Primofiore et al. (2020). In order to provide general results and some new insights with respect to the topography-related site effects, we iterate the numerical evaluation of ground motion for three different versions of the Arquata del Tronto visco-elastic sub-surface model, with growing complexity. The first version is characterized by a homogeneous structure, the second one presents near surface alterations (weathering layer and a restricted soil basin over a flat area), and the third one exhibits a complex internal structure and corresponds to the one used in Primofiore et al. (2020). The differences between the seismic responses evaluated for each model are then analyzed in three different frequency bands of engineering interest. We also analyze the possible correlation between the amplification and the topographic curvature, which was proposed as a proxy for topographic amplification by Maufroy et al. (2015).

Description of site response in terms of response function 3 × 3 matrix
In our characterization of the site response to a plane wave with vertical incidence we follow Paolucci (1999) and employ 3D impulse response functions evaluated from 3D physics-based numerical simulations of seismic wave propagation.
The 3D response function consists of a 3 × 3 matrix with components h ij (x, t). The time series h ij (x, t) describes the i-th spatial component of the seismic time series that is observed at the site x when the area is invested by a vertically incident impulsive plane wave polarized in the j-th direction with unitary amplitude. We can express the ground motion y(x,t) = [y 1 (x,t), y 2 (x,t), y 3 (x,t)] at the site x for any vertically incident seismic input where * denotes convolution in time. In virtue of the convolution theorem, the expression can be evaluated as a multiplication in the frequency domain where we denote with uppercase letters the Fourier transform of the corresponding timedomain quantities denoted in lowercase letters. (1) If an ideal elastic medium was considered, we could numerically evaluate the terms h ij (x, t) as the three component solutions of the wave field at the topographic surface for unitary impulsive plane wave input polarized in the three spatial directions (NS, EW, and Z) and with vertical incidence at the base of the model. In the present case however, a viscoelastic medium is considered and the intrinsic attenuation (and dispersion) would imply a response function which is dependent on the depth at which the plane wave is excited. In order to get rid of this incongruence, we evaluated the 3D response functions by making use of the reference site concept. In our approach, schematized in Fig. 1, we evaluate the solutions for polarized plane waves also in a "reference" spatial domain, characterized by a homogeneous medium (having the properties of the bedrock in the model under investigation) and flat topographic surface (with elevation equal to the average of the studied case). The depth at which the plane wave is introduced is irrelevant, as long as it is the same in the two cases. If we denote the frequency domain solutions in the reference medium and in the investigated medium with Y ij ref and Y ij respectively, (with the index i denoting the spatial component of the solution and the index j denoting the polarization of the excited plane wave) we can estimate the components of the response function of the model under investigation in frequency domain as the quotients Fig. 1 Schematic representation of the adopted approach for the evaluation of the 3 × 3 response function matrix described in Eq. (3). We denote with u 1 (t) the waveform of the polarized plane wave entering in the spatial domain from below with vertical incidence, with y 11 (t),y 21 (t) and y 31 (t) the three components of the corresponding solution in a station on the topographic surface of the heterogeneous viscoelastic model (on the left) and with y 11 ref (t) the only non-zero component of the solution on the flat surface of the homogeneous model (on the right). The scheme is repeated for the two orthogonal transversal polarizations and the longitudinal polarization of the plane wave This definition of the response functions also relaxes the request of the unitary and impulsive character of the plane wave entering the spatial domain. We remark that the present description of site response is valid if the medium remains in the linear range of induced deformations. In the present study we neglect possible non-linear effects.

Computation of response function from 3D model with SPECFEM3D
Fundamental factors for a successful numerical simulation consist in the availability of a plausible 3D sub-surface model of the spatial domain with an adequate geophysical characterization of the geological formations and in the accurate application of a numerical method for the solution of the elastodynamic wave equation.
In the present work, we constructed a digital 3D sub-surface model of the studied area by means of the commercial software Geomodeller3D, which interpolates the available geological data with methods based on the compliance of stratigraphic rules and on the interpretation of geologic surfaces as isopotential surfaces of a scalar potential field (Calcagno et al. 2008). We considered the geological and geophysical data from the surveys performed in the area during microzonation studies (ISPRA 2018;Puzzilli et al. 2019). Details on the 3D model we built are discussed in section Models.
Considering the need to solve accurately the wave propagation in a spatial domain bounded by an irregularly shaped free surface, we decided to perform the numerical simulations with the Spectral Element Method (SEM). This approach, first introduced in the study of acoustic and elastic waves by Seriani and Priolo (1994), combines the geometrical flexibility of Finite Elements Methods (FEM) with the accuracy of the spectral methods. In order to perform the simulations illustrated in the present work, we relied on the software SPECFEM3D Cartesian from source (Peter et al. 2011). The code accounts for the intrinsic attenuation properties of the viscoelastic medium by means of a series of standard linear solid elements (Savage et al. 2010) and avoids the boundary reflection artefacts thanks to the usage of the perfectly matched layer technique (Xie et al. 2014).

Reduction of the 3 × 3 response function to the plane wave amplification function
In seismic hazard studies the site effects are typically described as the ratio between the Fourier Amplitude Spectrum (FAS) |Y hor (x,f)| of the combined horizontal components of the ground motion at the observation site x and that of the seismic input |U hor (f)| where, having assumed that components 1 and 2 lay in the horizontal plane, we put . (4) and an analogous expression for |U hor (f)|. From Eq.
(2) appears that the amplification function A(x,f) defined in Eq. (4) necessarily varies with the considered seismic input, also because of cross-coupling effects (Paolucci 1999). Therefore, in our construction A(x,f) behaves as a random variable. We can expect a log-normal distribution for the spectral amplitude ratios A n (x,f) corresponding to different seismic inputs denoted by the index n, with geometric average A * (x, f ) and the multiplicative or geometric standard deviation (GSD) * (x, f ) defined as In this study we estimate the site average amplification and its multiplicative standard deviation by considering a set of three component earthquake ground motion records as seismic inputs U (n) (f) and by evaluating the correspondent theoretical site recordings Y (n) (x,f) by applying Eq. (2). Since we consider the numerical solutions for vertically emerging plane waves, we call this estimation of the site amplification as plane-wave amplification. The reduction of the 3 × 3 response function to the site amplification function can be therefore obtained by replacing in Eq. (6) the terms where H ij are 3D response functions as defined in Eq. (3) and where we omit the frequency and space dependence from the notation for the sake of simplicity. The GSD considered in the results therefore depicts the amplification function's uncertainty, that is due to the variable contribution of the input motion components to the output motion's horizontal components also through the cross-coupling effects (Paolucci 1999).
In the present analysis we considered as seismic inputs the ground motion records of 200 events with magnitude in the range from M = 2.8 to M = 4.8 that occurred during the 2016-17 Amatrice seismic sequence (see Figs. S1 and S2 in supplemental information) and that were collected at the reference site station MZ75 that was located 9 km north of Arquata del Tronto in the Montegallo municipality. Following Laurenzano et al. (2019), MZ75 represents the best choice as reference site for the Arquata del Tronto neighboring areas, because of its flat microtremor horizontal to vertical spectral ratio, its lithology which consists of the outcropping bedrock in the area (i.e. the arenaceous member of the Laga Flysch Formation, as described later in this article in the "Geomorphological and lithological characters" subsection of the "Study case" section), and because of its position in a valley where topographic effects should be minimal. The magnitude range and the relatively short distance of the events implied a high signal to noise ratio in the considered frequency band.

Horizontal polarization analysis from the 3 × 3 response function
Site effects on pronounced reliefs are often correlated with ground motion directionality (Burjánek et al. 2014;Massa et al. 2014) and a polarization analysis of the site response is therefore appropriate. We limit our polarization analysis to the horizontal components and following the scheme of the plane wave amplification analysis we define the plane wave directional amplification A * (x, f ) as the geometric average of the ratio of the FAS of the ground motion component in the direction φ measured clockwise from the 2nd axis (North) and the FAS of the seismic input in the same direction The amplitude spectra of the seismic motion |S(f)| in direction φ is evaluated from the complex spectra of the two horizontal components S 1 (f) and S 2 (f) as In particular, we set S 1 = U 1 (n) and S 2 = U 2 (n) when evaluating |U φ |. The polarization analysis consists in the computation of the directional amplification function for azimuth angles from 0° to 175° using a 5° step.
The distribution of the horizontal polarization in the site response is then analyzed in terms of the direction φ max of the maximum directional amplification and with the polarization level defined as where A min and A max are the minimum and the maximum directional amplifications at each point. P is zero when the amplification is isotropic and tends to 1 for increasing variability in the directional amplification.

Rotation rate analysis from the 3 × 3 response function
In consideration of the growing interest for the implications of the differential ground motion in seismic hazard (Guidotti et al. 2018), we attempted to analyze how topography and the interaction between topography and subsoil heterogeneities may induce torsional excitation (i.e., rotation around the vertical axis) at the Arquata del Tronto site. Given the ground velocity time series in a grid of points x at surface, we follow the notation suggested in Evans and International Working Group on Rotational Seismology (2009) and describe the rotation around the vertical axis (torsional motion) by means of the rotation rate In the case the ground motion is due to the passage of a plane wave, we can infer from Eq. (11) a proportionality between the rotational motion and the time derivative of the translational displacement (Newmark 1969). If we denote with Θ 3 (f ) the Fourier spectrum of the rotation rate around the vertical axis and with Ü (f ) the acceleration spectrum of transversally polarized horizontal wave motion (SH waves and Love waves), we can write where s(f) denotes the horizontal slowness, i.e. the inverse of the frequency-dependent apparent horizontal phase velocity of the propagating wave. The scaling factor between the rotational rate and the translational acceleration changes slowly and monotonically with the frequency (e.g. Trifunac and Todorovska 2001). The validity of Eq. (11) and of its equivalent in terms of peak ground rotation and peak horizontal velocity has been confirmed in field observations (e.g. Spudich and Fletcher 2008) and reproduced in numerical simulations (e.g. Stupazzini et al. 2009) in topography-free cases. Here we discuss instead the implications of Eq. (11) in the case of topographic features for a vertically emerging plane wave.
A vertically emerging plane wave in a homogeneous medium with horizontal flat surface is characterized by s = 0 and thus implies no rotation. On the other hand, heterogeneities and topographic features would imply local variations of s and consequently may induce torsional motion as a site effect. Simple geometrical considerations (see Fig. 2), allow us to deduce that in the ideal case of a homogeneous elastic medium with shear wave velocity V S , the topographic surface with slope β would imply for a vertically emerging plane SH wave an apparent horizontal slowness with value If we substitute Eq. (13) in Eq. (12) we can predict the local (i.e., dependent on the position x) scaling factor A θ between the Fourier amplitude spectra of the rotation rate and the acceleration of the seismic input as

Fig. 2
The relation between the horizontal slowness s and the shear wave velocity V S of a vertically emerging plane wave for a generic topographic surface with slope β (Eq. 13) can be deduced by considering the delay t BA between the wavefront's arrival times at two points A and B with different elevation An additional layer of constant thickness above the homogeneous bedrock (e. g. the weathering layer with shear wave velocity V S1 < V S ) implies a constant delay in the arrivals of the wavefront on the topographic surface in respect to the homogeneous case and therefore it does not change the apparent horizontal slowness of the SH wavefront on the topographic surface, which remains dependent only on the slope angle β and V S in the bedrock (Eq. 13).
On the other hand, the additional layer gives rise to Love waves which propagate in the direction parallel to the inclined topographic surface with a frequency dependent phase velocity c L (f), with values decreasing monotonically from V S in the bedrock at the low frequency limit, down to V S in the top layer, at the high frequency limit. In contrast to Eq. (13) the apparent horizontal slowness of Love waves on an inclined plane results Given that it is not trivial to predict the Love waves' amplitude from the amplitude of the incident wavefront that generates them on an irregular topographic surface, a strict prediction of the scaling factor between the rotation rate related to Love waves on a slope and the amplitude of seismic input is not practical. We can however expect that the scaling factor due to the contribution of Love waves to the torsional motion would behave as Given that c L (f) < V S and that tan β < (cos β) −1 it is indeed possible that rotation rates related to Love waves in the layered medium overcome the rotation rates related to SH-waves. In order to evaluate the effectiveness of Eq. (14) in the case of our 3D models, we apply Eq. (11) to compute the rotation rate time series ̇( n) 3 (t) from the horizontal components y 1 (n) (x,t) and y 2 (n) (x,t) of the ground velocity on the topographic surface, that are in turn estimated from the n-th seismic input using Eq. (1). Finally, we evaluate the plane wave torsional scaling factor A * as the geometric average of the ratio between the Fourier amplitude spectra of the rotation rate and of the horizontal acceleration of the seismic input and the corresponding multiplicative standard deviation. In fact, in analogy with the amplification function, the cross-coupling effects produce dispersion also in the scaling factor A θ and we can expect a log-normal distribution for A θ values estimated from different seismic inputs.

The case study
We used numerical modeling to estimate the site response in a dense grid of points distributed over a limited area centered on the location of the Arquata del Tronto village in Central Apennines (Italy).

Motivation
The extensive damage caused by the 2016 Central Italy seismic sequence in Arquata del Tronto village, set on a rocky asymmetric ridge and, on the contrary, the lighter damage observed in Borgo hamlet, set in a valley approximately 200 m north of Arquata, suggested that topography related modifications of earthquake ground motion took place in the area ( 2019) found that the spectral amplification at Borgo is characterized by a clear peak centered at 4.5 Hz, whereas the Arquata Castle presents broadband spectral amplification and frequencydependent polarization effects. Primofiore et al. (2020) performed 2D and 3D numerical simulations of seismic waves propagation using subsoil models of the Arquata del Tronto area that were deduced from the available geognostic, geophysical and geomechanical surveys. They found that 2D simulations, that are nowadays commonly used in microzonation practice (Vessia et al. 2021), cannot reproduce the seismic response of such asymmetric ridge, whereas 3D simulations that included the coupling effect of stratigraphy and topography, associated with the presence of a weathered portion of the rock mass, deposits, and the alternation of highly dipping rocky materials, provided results that are consistent with the observed heterogeneous damage distribution and with the seismic response characteristics at Borgo and Arquata Castle. Moreover, the snapshots on the topographic surface of the 3D seismic wave-field simulated in Primofiore et al. (2020) suggested that accidental eccentricity due to torsional ground motion may also be considered among topography-related site effects. Altogether, the framework that was set-up for 3D numerical simulations in Primofiore et al. (2020) permits to evidence effects on the seismic wavefield such as amplification, polarization and rotation at frequencies of engineering interest. Therefore we decided to exploit it to further deepen the understanding of some aspects of the topography's role in site effects.

Geomorphological and lithological characters
The geological substratum of the area is represented by the torbiditic Laga formation consisting of an alternation of arenaceous and pelitic-arenaceous lithofacies (Giallini et al. 2020;Amanti et al. 2020). The Laga formation is characterized by different stratified lithotypes, roughly 50° towards WNW dipping, and locally covered by fluvial deposits or arenaceous/calcareous debris (NE and SE area at the base of the Arquata crest). Instead, the lithotype beneath Borgo, in the valley below the hill, presents 10 to 30 m thickdeposits. In the description of the deep structure of the velocity model we consider two units (i.e. the pelitic and the arenaceous unit), whereas other two units are considered to describe the overall superficial 20 m thick alteration layer and the deposits at Borgo. Further details regarding the set-up of the 3D geological model are given in Primofiore et al. (2020). We illustrate the morphological and lithological characteristics of the area with the four maps in Fig. 3. A vertical cross section along the Arquata crest is shown in Fig. 4.
The elevation map evidences how the WNW-ESE oriented ridge of Arquata is transversely cut by saddles. The elevation stands in a range between 600 and 860 m a.s.l. In the geological units map we represent the distribution of the two deep geological units below the overall top alteration layer. In the same map we evidence the location of three sites of particular interest: Site A corresponds to the Arquata historical city center, Site B corresponds to the Borgo area in the valley and Site C corresponds to the Castle settlement on the crest. Site A and C are both placed on the Arquata ridge and lay above the same arenaceous formation, whereas site B, in the valley, mostly lays on the pelitic association with lower velocity (the gravel deposits that characterize site B are not shown in the map). The other two maps show the distribution of the topography's slope and curvature. If we locally approximate the elevation z(x,y) with the polynomial (Zevenbergen and Thorne 1987): The slope β and the curvature χ are defined (Moore et al. 1991) The curvature can be used as a general measure of the land convexity, in the sense that the positive values indicate the convex surfaces (ridges and peaks), and the negative values identify concave surfaces, i.e. valleys.
The steepest slopes of the ridge are found on the East and North sides of the promontory with a slope angle of 60 degree, whereas we can identify the area of high convexity with the main WNW-ESE ridge and its NS oriented secondary ridges. Comparing the geological units map and the curvature map, we observe that the NS oriented secondary ridges overlap the outcrop of the arenaceous formation.

Numerical models
We simulated the wave propagation in a spatial domain that represents a square area centered on the Arquata del Tronto hill, with 1 km long sides oriented along the geographical axes and that extends almost 1 km in depth.
(18) z(x, y) = Ax 2 y 2 + Bx 2 y + Cxy 2 + Dx 2 + Ey 2 + Fxy + Gx + Hy + I, Computer-aided 3D geological modeling allowed us to partition the domain in four main geological units (Table 1, AdT2 model), consistently with the geological map and a number of geological vertical sections (ISPRA 2018). The physical properties assigned to each geological formation were estimated from geophysical investigations conducted in the area in the aftermath of the 2016 earthquake (Puzzilli et al. 2019). We described the intrinsic attenuation of seismic waves with the shear quality factor Q μ , which was evaluated from V S with the widely used rule of thumb Q μ = 0.1 V S [m/s] (e. g. Paolucci et al. 2015) and with the bulk quality factor Q κ , which was in turn estimated as Q κ = 3.5 Q μ in accordance with the theory exposed by Morozov (2015).
In order to perform the numerical simulations we embedded the described 3D domain in a larger volume with 3 km long side square base aligned with the cardinal directions, extending from a depth of 500 m b.s.l. to the top surface which honors the 10 m resolution DEM Tinitaly (Tarquini et al. 2007). We discretized the spatial domain with an unstructured mesh, made of 518.400 hexahedral elements, using the SPECFEM3D internal mesher. We assigned to each element of the mesh the physical properties corresponding to the geological unit in its position. The elements' size is approximately 25 m near the top of the volume, where the low seismic velocities that characterize the alteration layer and alluvial deposits imply shorter wavelengths. On the other hand, we set an element size of 50 m in the lower part of the volume, where seismic velocities are higher. A 75 m thick Perfectly Matched Layer was applied on the lateral and bottom boundaries of the volume in order to exclude artefact reflections from the numerical solution. The whole numerical mesh is shown in Fig. 5a, whereas the part clipped on the central 1 km 2 study area is presented in Fig. 5b. The topography outside the central area was tapered towards the domain borders in order to exclude effects of topographic features not belonging to the central area.
To investigate the separate contribution of stratigraphic and topographic effects on the seismic ground motion, we carried out our analysis on three different geophysical models (Table 1 and Fig. 5): AdT0, AdT1 and AdT2. AdT0 presents a homogeneous structure and is therefore characterized by topographic effects alone. AdT1 includes a 25 m thick weathering layer below the surface and the alluvial deposits in the small basin at site B whereas AdT2 includes also the distinction between the arenaceous and pelitic associations in the deeper structure. The three models serve as canonical examples with different levels of complexity. Since it incorporates the available knowledge on the site geology, the model AdT2 can be indeed regarded as the most realistic, even though it was not strictly validated, as pointed out in Primofiore et al. (2020).
We modelled the vertically incident plane waves by placing a horizontal carpet of sources at 400 m b.s.l, with an explosive mechanism for the vertically polarized plane wave and West-East or North-South oriented double couple mechanism for the two independent cases of horizontally polarized plane waves, respectively. The homogeneous distribution of synchronized sources on a horizontal plane generates two horizontal wavefronts, one propagating downwards, which is absorbed in the bottom PML and the other upwards, which we use as the vertically incident plane wave required for the evaluation of the response function. The disturbance from the edges of the source surface is avoided by extending the source surface well outside the area of interest (that is also why we extend the spatial domain over a 3 km × 3 km spatial domain but focus only on the central 1 km × 1 km area) and by muting the solutions behind the theoretical arrival from the nearest edge of the source. With an analysis of the snapshots of the computed wavefield we verified that there is no superposition between the response to the incident plane wave and the disturbance from the edges at any receiver in the area of interest and that the PML works effectively. A set of three simulations is necessary to fully define the 3D response function. We Considering the values of the mechanical properties (Table 1), the dimensions of the mesh elements imply a time step as short as 0.0005 s to satisfy the stability criterion of the time marching step implemented in the SPECFEM3D solver and according to well established criteria in spectral elements usage (Seriani and Priolo 1994), accurate solutions up to the frequency of 9.6 Hz can be obtained. In our analysis we keep an additional margin of safety by rounding the maximum frequency to 8 Hz. On the other hand, given that the model does not contain any realistic feature that could meaningfully affect the seismic wave propagation at frequencies below 1 Hz, we neglected the solution at frequencies lower than that. Our case study was therefore focused on the seismic response in the threeoctave spanning frequency band 1-8 Hz.
Each simulation consisted therefore in 16,000 time iterations with a computational cost amounting to approx. 900 CPU hours on 225 cores on the Galileo supercomputer at CINECA.

Results and discussion
In order to identify the possible site effects that are due to surface morphology alone, we first computed the numerical response functions (see Eq. 3) for a homogeneous medium (model AdT0 in Table 1 and in Fig. 5) and evaluated over the topographic surface the geometric averages of the site response parameters defined in the Method section. Successive computations regarding inhomogeneous media (models AdT1 and AdT2 in Table 1 and in Fig. 5) allowed us to estimate the possible aggravations induced by the interaction between topography and stratigraphy. In order to evaluate the meaningfulness of the aggravations in the spectral amplification due to the changes in the model AdTx in respect to the reference model AdT0, (where AdTx stays for AdT1 or AdT2), we estimated the effect size, defined as with A ln and ln defined in Eq. (6). In consideration of the frequency dependence of the site response, we summarize the results with maps representing the arithmetic average value of the site response parameters functions over defined frequency bands. In particular we considered the arithmetic averages over the following frequency (period) bands: 1-2 Hz (0.5-1 s), 2-4 Hz (0.25-0.5 s) and 4-8 Hz (0.125-0.25 s).

Amplification
In the homogenous medium hypothesis (model AdT0), the topographic features at Arquata del Tronto generate amplification effects in all three considered frequency ranges, as can be observed from the maps shown in Fig. 6. As expected, the amplification at a given frequency is controlled by topographic features with characteristic sizes comparable with the associated wavelength (Geli et al. 1988). The average amplification computed in the 1-2 Hz band correlates well with larger topographic features such as the whole Arquata Because of flat topography site B is unaffected by amplification in the considered bands, whereas sites A and C appear most affected in the 2-4 Hz and 4-8 Hz bands, respectively. The amplification however just exceeds the value 2 in the latter case. The geometric standard deviation is relatively low in the valley and along the top of the ridge in a strip whose width progressively decreases with frequency and almost vanishes in the high-frequency band.
The qualitative correspondence between topographic features and the response function amplitudes suggests that typical topographic effects could be speedily estimated from local morphometric variables. A promising approach in this sense is the Frequency-Scaled Curvature (FSC) introduced by Maufroy et al. (2015), who found a correlation between the seismic ground-motion amplification estimated with the Median Reference Method MRM (Maufroy et al. 2012) and the curvature χ (defined in Eq. 16), with a wavelengthdependent smoothing applied on it. In analogy to Maufroy et al. (2015) we compared the plane wave amplification with the smoothed curvature for different values of the smoothing length. We used the SAGA algorithm (Conrad et al. 2015) to evaluate the curvature χ at the surface points of the spatial domain and applied a n × n smoothing kernel, where n is an odd number which defines the smoothing length as L s = 2n Δx where Δx is the spatial increment (25 m in our case).
Our comparison consists in the evaluation of the linear correlation coefficient r between the average value of the natural logarithm of the amplification function in the three Fig. 7 (Top) correlation coefficients between logarithmic plane-wave amplification ⟨A ln ⟩ f 1−f 2 and smoothed curvature χ in three frequency bands for different values of the smoothing length Ls; bottom) weightedleast square regression analysis between ⟨A ln ⟩ f 1−f 2 and smoothed χ with smoothing length Ls corresponding to maximum correlation for each frequency band frequency bands ⟨A ln ⟩ f 1−f 2 and the smoothed curvature χ (Ls) , where we express the smoothing length L s in terms of the S-wavelength λ c , defined as the geometric mean of the wavelengths corresponding to the band extremes. From the plots in the upper part of Fig. 7 it appears that the maximum correlation corresponds to a smoothing length L s between 0.5 λ c and 0.8 λ c . This result is consistent with Maufroy et al. (2015), who established a maximum correlation between MRM amplification and smoothed curvature for a smoothing length of about half the wavelength. In the lower part of Fig. 7 we represent the correlation between ⟨A ln ⟩ f 1−f 2 and χ (Ls) with L s corresponding to the maximum correlation case with the best fit line of the form The coefficients a and b in Eq. (21) were evaluated for each frequency band and are reported in Table 2.
Given the simplicity of the AdT0 case (homogeneous medium), the increased dispersion level in the distribution of the amplification in relation with the smoothed topographic curvature in the 2-4 Hz (R 2 = 0.57) and 4-8 Hz (R 2 = 0.54) frequency bands can be only associated with issues regarding the ratio between the spatial sampling step and the wavelength. We have to consider that the evaluation of the curvature is based on the local polynomial approximation of the topographic surface as expressed in Eq. (18), whereas our numerical computation of the wavefield is based on the local approximation of the topographic surface with inclined flat quadrangles. The mismatch between the two approximated surfaces (and between each of them and the targeted real world surface), can be made negligible by choosing a sufficiently small spatial sampling step length in respect to the considered wavelength. Maufroy et al. (2015) considered as much as 70 sampling steps per wavelength in order to evidence the correlation between curvature and topographic amplification. Here we find that 20 sampling steps per wavelength (the 1-2 Hz case) allow for a coefficient of determination as high as R 2 = 0.8, but for smaller wavelengths (the 2-4 Hz and 4-8 Hz cases imply 10 and 5 sampling steps per wavelength, respectively) the scattering in the amplification vs. curvature plots appears significant, even though the correlation is still noticeable (Fig. 7, bottom).
In Fig. 8 we can see how the horizontal amplification functions in the three sites A, B and C change because of the introduction of complexities in the model structure. The most remarkable changes concern site B, which is characterized by flat topography. The introduction of soil basin increases the amplification at site B to higher values than the topography-related amplification at sites A and C. It is also evident from Fig. 8, that for frequencies above 2 Hz, the impact of model AdT2 is rather limited in comparison to model AdT1, a fact that can be explained with the relatively low impedance ratio of 1.14 between the arenaceous and pelitic associations in comparison of the impedance ratio of 1.56 assumed between the arenaceous bedrock and the weathered layer. In the maps on the left side of Fig. 9 we represent the spatial distribution of ⟨A * ⟩ f 1 −f 2 for the model AdT1, and on the right side the corresponding effect size in respect to the model AdT0.
The AdT1 model differs from AdT0 in a 25 m thick weathering layer just below the surface over all the area and a 10 to 30 m thick layer of gravel deposits in the alluvial basin in the restricted area of Borgo (site B).
The weathering layer leads to a general increase of the amplification, with modulations that follow the topographic features. In the low frequency band (f < 2 Hz) the increase in the amplification is not significant at any point, whereas in the intermediate and high frequency bands the amplification gains significance at some spots, but not on the flanks of the hills. The intermediate and high frequency amplification associated with the introduction of the soil layers in the alluvial basin at site B exceeds any amplification related to the interaction between topography and the weathering layer. The absence of amplification aggravation in the low frequency band can be easily explained if we consider that we assigned Vs = 450 m/s to the soil in the basin with a maximum depth of 30 m. In fact, in the described case the expression for the fundamental frequency of one-layer structures (Haskell 1960) predicts stratigraphic amplification effects for frequencies above 3.75 Hz, whereas if we combine the basin and the weathering layers in a single layer with average Vs, the resulting stratigraphic amplification does not affect frequencies below 2.5 Hz.
In Fig. 10 we represent an analogous comparison between the results concerning model AdT2 and model AdT0. The increase in the amplification that is due to the introduction of the deep structure (model AdT2) presents a pattern that is quite similar to the model AdT1, but with an overall effect size which is only slightly larger than in the AdT1 case.
In order to emphasize the differences in the effects of the models AdT2 and AdT1 compared to the model AdT0, we plot in Fig. 11 a histogram of the effect size values in the points on the map in the three considered frequency bands. The models AdT2 and AdT1 present differences in the amplification in the lowest frequency bands (1-2 Hz and 2-4 Hz), whereas in the high frequency band (4-8 Hz) the two models are essentially equivalent.  The broader distribution of the effect size values for the high frequency band in Fig. 11, can be explained with the sensitivity of the seismic response at high frequencies to small scale topographical and geological features.
With the introduction of inhomogeneity in the model we found a frequency dependent drop in the correlation between the smoothed curvature and average amplification. In Fig. 12 we plot the correlation scores for all the three investigated models in the three frequency bands. We note that in the 1-2 Hz band the heterogeneity in the model has almost no effect on the correlation between the amplification and the smoothed curvature. On the other hand, the variations in the subsurface model reduce substantially the correlation in the other two frequency bands, with the 4-8 Hz band being the most sensible to the subsurface structure effects, but without significant distinction between the two heterogeneous models. The observed additional disruption of the correlation (it adds to the one discussed for the homogeneous case), can be in part attributed to the stratigraphic amplification, which comes into play in the 2-4 Hz and 4-8 Hz bands in the considered case, and in part to the change of the effective wavelength caused by Fig. 11 Histograms of the amplification effect size values over the grid points sampling the studied area for the frequency bands: 1-2 Hz, 2-4 Hz and 4-8 Hz. The median value is indicated in order to assess the overall differences between the changes in the amplification functions due to the variations in the sub-surface models Fig. 12 Correlation coefficients between averaged horizontal plane-wave amplification and smoothed curvature obtained in the three frequency bands of interest for homogeneous model AdT0 (black) and inhomogeneous models AdT1 (dark grey) and AdT2 (light grey) the addition of a low velocity layer on the topographic surface, which would require an adaptation of the smoothing length in respect to the homogeneous medium case. As described in Wang et al. (2018), the most convenient smoothing length is a function of the layer thickness. For the 1-2 Hz band the weathering layer thickness is of the order of only 5% of the considered wavelengths and therefore its effect on the correlation between amplification and curvature appears minimal.

Horizontal polarization
The distribution of the horizontal polarization in the site response for the model AdT0 is summarized on the maps in Fig. 13, with the headless arrows indicating the direction φ max of the maximum directional amplification averaged over the three frequency bands and the arrows' length indicating the polarization level as defined in Eq. (10).
We can observe that the direction and level of polarization change significantly in the three considered frequency bands. If we compare the maps in Fig. 13 with the amplification maps in Fig. 6, the amplification appears to be generally uncorrelated with the polarization level: in many points the polarization of the response corresponds to directional de-amplification rather than amplification. This is the case of the polarization level found along the northern and the southern slope of the Mt. Arquata ridge in the 1-2 Hz and 2-4 Hz bands (maps on the left and in the center of Fig. 13). In the 2-4 Hz band (central map in Fig. 13), the polarization level appears to be larger in the portion of the ridge's flanks that is closer to the top, but actually no clear relation between the polarization direction and the topographic features can be observed at these points, where polarization does not imply amplification. On the other hand, the points that are characterized by both polarization and amplification are essentially located on the primary and (in the 4-8 Hz band illustrated in the bottom-left map in Fig. 6 and in the map on the right in Fig. 13), secondary ridges. In this case the polarization direction is typically perpendicular to the crest orientation. We can therefore conclude that even if the topography-induced amplification is usually associated with polarization, the polarization itself does not necessarily implies topography-induced amplification. In the same manner, the absence of polarization does not necessarily imply the absence of topographic amplification, as can be seen in the central points at site C in the 4-8 Hz frequency band, map on the right in Fig. 13). The lack of correlation between amplification and the polarization level is also clear from the plots in Fig. 14, showing the distribution of the amplification values in function of the polarization level for the homogeneous model.
The introduction of the weathering layer and deeper geological units (see Fig. S3 in supplemental information), does not bring noticeable changes in the polarization pattern in the 1-2 Hz and 2-4 Hz bands (top and center maps in Fig. S3), whereas in the 4-8 Hz band (bottom maps in Fig. S3), we can observe a reduction of the polarization level at site A on the edge of the ridge and the appearance of polarization in the seismic response at the edges of the nearly circular soil basin near site B, with direction tangential to the basin edge. The heterogeneous structure may therefore attenuate the topography-related polarization effects (without reducing the amplification) on one side and can produce topography-unrelated amplification and polarization, on the other. As can be seen from Fig. 15, the correlation coefficient between the amplification and polarization level for the heterogeneous models remains as low as it is for the homogeneous model. Fig. 13 Distribution of the horizontal polarization computed for the homogeneous model in the three frequency bands 1-2 Hz, 2-4 Hz and 4-8 Hz. Grey headless arrows in the maps indicate the maximum polarization direction according to the definition provided in the "Horizontal polarization analysis from the 3 × 3 response function" section. The polarization level, as defined in the same section, is indicated by the arrow length 1 3

Torsional motions
In Fig. 16 we map the peak rotation rate Θ max around the vertical axis evaluated from the seismic response of the three models for the seismic input corresponding to a M 4.8 event at a distance of 24 km. The maps demonstrate the possibility of topography induced torsional motions for vertically incident plane waves. The order of magnitude of the induced  Values of the linear correlation coefficient between the polarization level and amplification for the three subsoil models in the three considered frequency bands Fig. 16 Maps of the peak rotation rate around the vertical axis evaluated in the three models for the specific case of the seismic input corresponding to a M 4.8 event at a distance of 24 km torsional motions (~ 0.1 mrad/s) is compatible with the ones observed or predicted in equivalent cases in bibliography (e. g. Spudich and Fletcher 2008). The amplification effect of the introduction of heterogeneous subsoil on the topography-induced torsional motion is also evident.
In the following we analyze the possibility of induced torsional motion in the study area by means of the parameter A * defined in Eq. (14), which describes the scaling factor between the site-induced angular velocity around the horizontal axis (i.e. torsional motion) and the translational acceleration of the vertically incident plane wave. In Fig. 17 we plot, for the three models, the maps corresponding to the A * averaged in the 4-8 Hz frequency band. From the map corresponding to the homogeneous medium (model AdT0), it appears that the highest values of the induced rotational motion are associated with the points on topographic slopes, as predicted by the conceptual model that lead us to Eq. (13). The variability in ⟨A * ⟩ 4−8Hz is described by a multiplicative standard deviation which ranges between a factor 1.5 and 2. Maps regarding the frequency ranges 1-2 Hz and 2-4 Hz are available in the supplemental information as Figs. S3 and S4. From the maps regarding model AdT1 and model AdT2 in the frequency range 4-8 Hz (Fig. 17), it is clear that the introduction of the weathering layer enhances the rotational motion on the slopes whereas the introduction of the basin structure produces rotations near the basin edges. These basin-induced rotations are of the same order of magnitude as those induced by the concurrence of topography and weathering. We can also notice that since they are characterized by low slope angles, the sites A and C are not affected by the appearance of torsional components in the ground motion.
To verify quantitatively if Eq. (13) could be used to estimate the rotational scaling factor ⟨A * ⟩ 4−8Hz computed from the three models, we apply the weighted least squares approach to find the parameters a and b in the linear fit We consider the slope β distribution evaluated from the topography with the SAGA algorithm (Conrad et al. 2015), and apply a smoothing with length 125 m, corresponding to the wavelength of the incident waves in the AdT0 model at the maximum considered frequency (8 Hz). In analogy with the curvature-amplification case, the chosen smoothing length maximizes the linear correlation coefficient between the logarithms of the slope's tangent and the rotational scaling factor. We neglect the points with slope β < 15° in order to exclude from the analysis the rotation induced in the basin area near site B. The linear fit is shown in Fig. 18 and the parameters a and b corresponding to the three models are tabulated in Table 3 along with the V S value interpreted from the b parameter.   (13). The AdT1 and AdT2 models instead provide ⟨lnA * ⟩ 4−8Hz values that are less correlated with the topographic slope. As predicted in the subsection "Rotation rate analysis from the 3 × 3 response function" of the"Method" section, a more sophisticated predictive model than the one provided in Eq. (13) is required in the cases where the topography effect combines with stratigraphic effects.
In Fig. 19 we analyze the differences between the induced torsional motion effect size corresponding to the substitution of the AdT0 model with models AdT1 and AdT2, respectively. Only in the 4-8 Hz frequency band the median value of the effect size approaches unity, but in this frequency band the differences in the effect size of the two models are minimal. As in the amplification case, the 4-8 Hz band shows to be the most sensible to the subsurface structure effects, without significant distinction between the two heterogeneous models. Even though the two heterogeneous models bring different levels of size effects in the other two frequency bands, these levels are below unity and could be considered less significant. The homogeneous model assumption is therefore sufficient to evaluate induced rotations in the frequency band up to 4 Hz.

Conclusions
In the present work we aimed to provide some general results and new insights with respect to the topography-related site effects, by further exploiting the physics-based modelling of 3D seismic wave propagation in the Arquata del Tronto area, that had been set up in Primofiore et al. (2020). In particular we analyzed the separated and combined effects of irregular topography and subsurface heterogeneity on the seismic ground motion in terms of horizontal component spectral amplification, including the aspects that concern polarization and rotation. To facilitate the discussion, we averaged the results over the frequency (period) bands of engineering interest: 1-2 Hz (0.5-1 s or low frequency), 2-4 Hz (0.25-0.5 s or intermediate frequency) and 4-8 Hz (0.125-0.25 s or high frequency). By comparing the numerical results based on three different geophysical models of the area, Fig. 19 Histograms of the induced torsional motion A * effect size values over the grid points sampling the studied area for the frequency bands: 1-2 Hz, 2-4 Hz and 4-8 Hz. The median value is indicated in order to assess the overall differences between the changes in the amplification functions due to the variations in the sub-surface models one homogeneous, one with a weathering layer and one with a complex 3D structure, we were able to distinguish typical topographic effects from combined effects of topography and stratigraphy. We evaluated the intrinsic variability in the 3D seismic response caused by cross-coupling of different components of the ground motion by means of the sample standard deviation over a seismic input population composed of 200 three-component ground motion records collected in a station on seismic bedrock in the same geographical area. We used the evaluated standard deviation to estimate the effect size of the adoption of models with a heterogeneous subsoil in place of a homogeneous one.
The obtained results highlight the frequency dependence of topography-related amplification, which is controlled by topographic features with characteristic sizes comparable with the associated wavelength. Following an approach introduced by Maufroy et al. (2015), we found for the homogeneous model a good correlation between the amplification and the smoothed topography's curvature. In the two higher frequency bands we noted however a progressive scattering between the curvature and the amplification. This scattering can be attributed to the increasing mismatches between the effective topographic surface and the polynomial surface used to evaluate the curvature. In the analyzed case it appears that a sampling density of the order of 20 points per wavelength is required to evaluate curvature values for reliable estimation of the topographic amplification.
The maximum amplification values in a homogeneous model are however moderate (they do not exceed a factor of 2.5-see Fig. 6), but the introduction of the heterogeneous structure causes an aggravation, with modulations that follow the topographic and geological features. In the intermediate and high frequency bands, the aggravation is significant at some spots, in particular in the flat area occupied by the small alluvial basin, where the amplification equals and even exceeds the topography-related amplification at topographic highs (see Figs. 9 and 10). Because of the contribution of the stratigraphic effects to the amplification, the correlation between amplification and curvature further weakens for frequencies above 2 Hz, when heterogeneous subsoil is considered (see Fig. 12).
The directional analysis of the amplification in the homogeneous model shows frequency dependence, as well. A visual analysis of the polarization level and angle distribution over the area (see Fig. 13), reveals that on the topographic reliefs, where the polarization is associated with convexity-related amplification, the polarization tends to be orthogonal to the topography major axis in the low frequency range, whereas in the high frequency range, the polarization changes and tends to be orthogonal to secondary ridges axes. The same analysis reveals that many points on the map are characterized by a polarized response that does not imply amplification. Actually we found that the polarization level is generally uncorrelated with the amplification (see Figs. 14 and 15). We conclude that the polarization itself does not necessarily indicate amplification, since it could be associated with concavity-related de-amplification. A heterogeneous substratum seems to either attenuate the topography-related polarization effects or produce topography-unrelated amplification and polarization (Fig. S3 in the supplemental material).
Finally, we analyzed the possibility of topography-induced torsional motions in terms of the scaling factor between the induced rotation rate and the seismic input acceleration. From simple geometrical considerations we expected a linear proportionality between the scaling factor and the tangent of the topographic slope, which was confirmed through a linear regression of the results concerning the homogeneous model case. On the other hand the introduction of the weathering layer clearly implies an aggravation of the induced torsional motion, especially in the high-frequency band (Fig. 17, bottom-left map), which contains the layer's resonance frequencies. The torsional motion induced by the concurrence of slope and weathering appears to be of the same order of magnitude as that induced by the subsoil heterogeneity on the basin edge. This finding suggests that induced torsional motion could concur in the aggravation of earthquake effects on rocky hills like the one on which Arquata del Tronto rests.
For both the translational and torsional components, the most evident aggravation of topography-related effects in respect of the homogeneous model is given by the introduction of the weathering layer in the frequency band 4-8 Hz (rightmost panel in Fig. 11 and in Fig. 19). On the other hand, the median aggravation associated with the 3D structure exceeds the median aggravation associated with the weathering layer in the band 1-4 Hz, but it appears relatively less significant if measured in respect of the variability of the 3D seismic response related to the components' cross-coupling.
For simplicity, only the case of vertically incident plane waves was considered in the present work. The issue of possible modifications in the seismic response for oblique wavefronts in function of the back-azimuth and the incidence angle requires further investigation.