Introduction
Topography effects have been extensively recognised during strong ground motions. Site response is significantly affected around topographic irregularities and consequently a larger concentration of damage has been reported near topographic reliefs.
A number of previous numerical studies on topographic effects focus on a parametric assessment of the different geometric variables that are shown to affect the topographic amplification factors and mostly consider homogenous soil stiffness profiles, which are not completely representative of real site conditions.
This article aims to investigate the effects of canyon topography on ground motion for different soil stiffness profiles over a rigid bedrock using the finite element method. The soil stiffness profiles are selected to have similar oscillation characteristics to the uniform stiffness profile. The methodology for the selection procedure and the associated assumptions are detailed below. The numerical results herein focus on the variation of the maximum normalised acceleration with the examined period of input motion for all the studied profiles.
Background review
Topography effects are closely related to the presence of strong topographic reliefs (e.g. slopes, hills and canyons), complex subsurface geometries (e.g. alluvial valleys) and lateral geological irregularities (e.g. ancient faults). There is sufficient observational evidence from destructive earthquakes to conclude that concentration of damage occurs where steep slopes or complex topography is noticeable. Also, buildings and structures located on the top of hills and canyons, have suffered more intense damage than others located at the base of the same formations. Some of these observations refer to the San Fernando earthquake in 1971 (^{Boore, 1972}), the Chile earthquake in 1985 (^{Celebi, 1987}) and the Northridge earthquake in 1994 (^{Bouchon and Barker, 1996}). More recently, ^{Rathje et al. (2011)} observed an important damage concentration at the ridgetops and steeper slopes after the Haiti Earthquake in 2010. ^{Hancox and Perrin (2011)} indicated similar damage concentration at cliff tops and low ridge crests during the 2010 – 2011 Canterbury earthquake sequence.
Numerical studies on topographic effects are mainly performed using the Finite Element Method FEM, the Finite Difference Method FDM, the Spectral Element Method SEM and the Boundary Element Method BEM, among other numerical techniques (^{Tripe et al., 2013}). However, most of the numerical studies have not been able to reproduce the high amplification amplitude values recorded in the field at various case studies (^{Pedersen et al., 1994}; ^{Geli et al., 1988}). These differences have been attributed to factors that have not been considered in the numerical models, such as the 3D shape of the topographic features, the presence of adjacent topography and the effect of soil layer amplification (^{Tripe et al., 2013}; ^{Ashford et al., 1997}; ^{Bouckovalas and Papadimitriou, 2005}).
Most of the previous numerical studies are parametric, focusing on the type of incident wave (^{Assimaki and Gazetas, 2004}); frequency content of input motion (^{Bouckovalas and Papadimitriou, 2005}); slope inclination (^{Ashford et al., 1997}, ^{Bouckovalas and Papadimitriou, 2005}; ^{Tripe et al., 2013}); effect of cycles of input motion and soil damping (^{Bouckovalas and Papadimitriou, 2005}); soil stratigraphy (^{Assimaki and Gazetas, 2004}) and depth of bedrock (^{Tripe et al., 2013}), among other parameters.
^{Ashford et al. (1997)} firstly proposed that topographic amplification could be examined independently to the soil layer amplification. However, more recent studies (^{Tripe et al., 2013}; ^{Assimaki and Gazetas, 2004}; ^{Rizzitano et al., 2014}; ^{Assimaki et al., 2005}; ^{Assimaki and Jeong, 2013}) showed that topographic amplification is dependent on soil layer amplification and these two phenomena cannot be examined separately. These studies also present a better representation of the recorded ground motion on site as more realistic soil conditions are considered in the more recent numerical studies. ^{Rizzitano et al. (2014)} showed that a smaller impedance ratio between the soil layer and bedrock results in larger interaction of soil stratigraphy and topography. This is also shown in ^{Assimaki and Jeong (2013)} and was first highlighted by ^{Tripe et al. (2013)} who considered a rigid bedrock assumption in the numerical investigation of topographic effects on soil slopes. The importance of the simulation of the soil stratigraphy was also highlighted by ^{Assimaki et al. (2004}, ^{2005a}) during the investigation of the topographic effects associated with the Kifisos river canyon banks in the Athens 1999 earthquake. The effect of a double soil layer presence in a halfspace has recently been considered by ^{Assimaki and Mohammadi (2017)} with the aid of a numerical parametric investigation on topographic effects for hill geometries. The underlying soil conditions are shown to be important for the topographic amplification variation and its magnitude. These factors are further investigated herein for canyon cases within soil profiles founded on a rigid bedrock.
Methodology
Model description
Twodimensional 2D plane strain time domain finite element FE analyses were performed, considering a canyon geometry in a soil layer deposit over rigid bedrock. The geometry of the FE model is presented in Figure 1. The geometrical parameters for the current set of numerical analyses are kept constant and equal to: height of the canyon slope H = 50 m, slope angle i = 75^{o}, distance to bedrock z = 125 m, cresttocrest distance L_{ctc} = 280 m and distance of the slope crest to the lateral boundary of the mesh L = 560 m. The input motions correspond to vertical propagating inplane shear waves (SV waves) the soil is modelled as a linear viscoelastic material with varying stiffness with depth.
The study was carried out using the Imperial College Finite Element Program ICFEP (^{Potts and Zdravković, 1999}), employing the generalisedα time integration method. In terms of boundary conditions BC, for the static part of the analysis the horizontal and vertical displacements are equal to zero at the bottom boundary of the mesh, while the horizontal displacements are zero on the lateral boundaries. The same boundary conditions to the static part of the analysis are imposed on the bottom boundary during the dynamic part of the analysis due to the rigid bedrock assumption, while normal and tangential dashpots are used on the lateral boundaries.
The spatial discretisation of the mesh follows the recommendations of ^{Kuhlemeyer and Lysmer (1973)} where the element dimension, Δl, was considered as Δl ≤ λ/10, where λ corresponds to the wavelength related to the highest frequency of the considered input motion, i.e. for T_{p} = 0.1 s and V_{s} = 500 m/s, λ = 50 m thus an element size of Δl = 5 m was used.
The domain reduction method DRM is used to ensure that freefield conditions are reached at the lateral boundaries of the FE mesh and simultaneously to reduce the computational cost of the analyses. The DRM is a twostep procedure that intends to reduce the domain of the problem by introducing changes in the governing variables. The DRM was first developed by ^{Bielak et al. (2003)} for seismological purposes and it has been extended and implemented in ICFEP for dynamic coupled consolidation problems (^{Kontoe et al., 2009}).
All the analyses have been performed using wavelets of a harmonic nature with predominant period T_{p} as input motions. This is a harmonic wavelet which is modulated by the ^{Saragoni and Hart (1973)} temporal filter and corresponds to the Chang's time history used in previous studies (^{Bouckovalas and Papadimitriou, 2005}). The expression used for the acceleration time history is given by equation (1)).
where α, β and γ refer to the constants controlling the shape and amplitude of the wavelet, T_{p} corresponds to the predominant period and t is the time. The values of each constant were varied for each considered T_{p} value to reach a maximum amplitude of unity. Additionally, the number of cycles of the input motion is maintained constant and equal to 12 for all the studied input motions. A plot of the acceleration time history is presented in Figure 2 for T_{p} = 1 s and constant values of α = 2.0, β = 1.5 and γ = 5.0. The considered input motions refer to an input motion period range T_{p} from 0.1 s to 2.0 s.
The Step I of the DRM comprises of a 1D soil column analysis with a soil thickness equal to the crest height of the 2D canyon model (i.e. equal to the depth to bedrock z value). This analysis represents the freeheld response and takes into account the soil layering effects. The results of this analysis are imposed on the Γ line of the 2D numerical model (see Figure 1) during the Step II of the numerical process. The model of Step II incorporates the local effects that refer to the presence of the topographic irregularity in this case.
1D analyses were performed for each predominant period T_{p} of the input motion and considered the following steps. The 1D model is generated taking into account the mesh recommendations (^{Kuhlemeyer and Lysmer, 1973}). The motion is applied in the horizontal direction along the base of the model, while tied degrees of freedom TDOF were applied between the nodes at the same elevation along the lateral boundaries. The dynamic analysis for the 1D crest model is used to generate the input for the Step II of the DRM, but also to obtain the freefield crest acceleration response. The freefield response is then used to normalise the horizontal and vertical 2D response obtained from the Step II analysis in terms of accelerations. Similarly, dynamic analyses for a 1D toe model (i.e. reduced depth to bedrock corresponding to the toe conditions) are performed in order to obtain the freefield acceleration at the toe.
2D analyses were also performed for each predominant period T_{p} of the input motion, considering the following steps. Firstly, 2D static analysis is carried out modelling the geostatic stress field for level ground and then part of the material is removed to simulate the canyon geometry in a single stage of excavation. Following this, the 2D dynamic analysis is performed, using as input the results from the DRM Step I 1D model. For each 2D analysis, acceleration time histories were recorded along different points on the surface of the canyon geometry.
Topographic amplification is defined as the ratio of the resulting motion close to the irregular feature (topo) to the motion in the freefield (ff). The motion at the canyon surface is characterised by both horizontal and vertical components. The vertical component is of a parasitic nature as it results from the wave scattering (refraction and reflection) at the irregular ground surface. With the aim of decoupling topographic effects from soil layering amplification, a normalisation process is used. The topographic factors in the horizontal and vertical direction result from the normalisation of the horizontal and parasitic vertical components of ground motion to the horizontal component of the freefield ground motion in the case of SV wave incidence. This is because there is no vertical component of the ground motion in the freefield as the incident waves are simply reflected at the flat ground surface with the absence of wave scattering.
The examined soil profiles
Three soil profiles are examined herein, one with a linear stiffness variation with depth (Linear Constant), a trilinear stiffness profile representing a case of a parabolic stiffness variation with depth (Toro V_{s}) and a stepped stiffness profile representing two soil layers with a strong stiffness contrast within the soil layer above the rigid bedrock (Stepped V_{s}). These cases were chosen because they correspond to more realistic spatial variations of soil stiffness encountered in nature. The Toro V_{s} profile is based on the Toro model (^{Toro, 1995}), which corresponds to a shear wave velocity profile generated by statistical approaches. In the present study, a lower bound profile was used, which corresponds to a median profile estimated with the Toro process minus one standard deviation, approximated by three linear parts.
All the examined profiles are derived so that they have a comparable first mode of vibration to a reference uniform profile with shear wave velocity of V_{s}= 500 m/s. This is to ensure that the soil layer amplification is similar for all the examined cases of stratigraphy, as the fundamental soil layer modes of vibration are seen to affect the frequency content of the topographic amplification (^{Tripe et al., 2013}; ^{Skiada et al., 2017}, ^{2017a}).
Two weighted average values of shear wave velocity are used herein: according to the ^{Dobry et al. (1976)} and as proposed in the ^{IBC (2003)}, ^{Eurocode 8 (1998)} and ^{ASCE 702 (2005)}. Both definitions are expressed in equation (3)).
Both average values of shear wave velocity are compared to the uniform stiffness profile, aiming to have a similar value. The considered stiffness variation with depth and the main characteristics of each profile are presented in Table 1 and Figure 3 in comparison to the uniform profile.
Parameter  Uniform  Linear constant  Toro V_{s}  Stepped V_{s}  

Shear modulus, G, MPa  500  100 + 6.4z  Varying in depth (see Figure 3)  125  z < 65 m  
1125  z > 65m  
Bulk modulus, K, MPa  1333  266.6 + 17.06z  Varying in depth (see Figure 3)  333.3  z < 65 m  
3000  z > 65m  

500  493  529  500  250  z < 62.5 m 
750  z > 62.5 m  
V_{s125}, m/s  500  459  479  500  
Fundamental period T, s  1  0.906  0.855  1.085 
z: depth from soil surface level (crest level) in meters
The soil mass density is equal to ρ = 2.0 Mg/m^{3} and the horizontal coefficient of earth pressure K_{0} = 1.0 in all the performed analyses. For the dynamic analyses, the target damping ratio, ξ = 5%, is achieved by varying the Rayleigh damping parameters according to the input motion period and the fundamental period of the canyon.
The transfer functions F_{2}(ω) resulting from the considered profiles are plotted in Figure 4 in comparison to the transfer function of the uniform stiffness profile. The transfer functions are calculated using the software STRATA (^{Kottke and Rathje, 2013}). The transfer functions of the varying soil stiffness profiles are characterised by higher amplitudes of F_{2}(ω) and shifts in the position of the natural periods compared to the uniform stiffness profile (observed at ~ 1 s, ~ 3 s, ~ 5 s, etc). The location of these shifts and the resulting amplitude values were confirmed theoretically using the expressions proposed by ^{Ambraseys (1959)} and ^{Dobry et al. (1976)}.
Results
Figure 5 shows the horizontal A_{hmax} and vertical A_{vmax} topographic amplification factors for the Uniform and Linear Constant profiles along with their respective Transfer Function F_{2}(ω) at the crest of the canyon x_{cr} = 0 m in Figure 1). The obtained results confirm that the topographic amplification maximum amplitude A_{hmax} is located between the first and the second natural period of the 1D profile. This effect could be explained as follows: the maximum normalised response occurs when the Transfer Function F_{2}(ω) reaches its minimum value, due to the normalisation process by the 1D soil column (^{Skiada et al., 2017}).
The horizontal A_{hmax} and parasitic vertical A_{vmax} topographic amplification at the crest of the canyon (x_{cr} = 0 m in Figure 1) with input wavelet period T_{p} are presented in Figure 6. It can be observed that the amplitude of both the horizontal and the parasitic vertical topographic amplification differs between the several stiffness profiles. Linear Constant and Toro V_{s} profiles are characterised by a peak of A_{hmax} = 2.5 at T_{p} = 0.5 s in the horizontal direction and on overall similar response for a wide frequency range. The similarity of the topographic response in these cases is expected because of the considered profile transfer functions which are similar and the stiffness variation which is mostly different in the first few meters below the ground surface. This difference mainly affects the higher modes of vibration of the 1D soil columns (see the larger difference of the examined profiles to the uniform one at the second and third mode in Figure 4). The difference of the topographic amplification is then expected to be seen in the higher frequency modes (i.e. in the smaller input motion periods T_{p} area). Overall, the horizontal topographic amplification is similar for both the Linear Constant and the Toro V_{s} profiles as well as the Uniform profile, with some differences in the high frequencies range.
The Stepped V_{s} profile is, however, characterised by a secondary peak at T_{p} = 0.16 s additional to the peak position at T_{p}= 0.5 s with A_{hmax} = 1.75. The different peaks of this response compared to the uniform profile are expected due to the presence of a double layer formation in this case. The top layer is also three times softer than the bottom one (impendence ratio = 3) in this case, so it is characterised by a more complex wave scattering mechanism compared to the uniform stiffness profile.
The parasitic vertical response of the Uniform profile is characterised by a maximum located at T_{p} = 0.5 s with magnitude value over A_{vmax} = 2.0. The Linear Constant and Toro V_{s} profiles also maximise at T_{p} = 0.5 s. However, a secondary peak at T_{p} = 0.160.2 s of similar amplitude to the peak at T_{p} = 0.5 s is a characteristic of these profiles. This highlights the largest effect of the change of the stiffness profile on the vertical response. In this case, the amplification of the vertical response is also comparable to the horizontal one. A similar response is observed for the stepped profile, as it is characterised by two peaks at 0.16 and 0.5 s, however with smaller amplitude in comparison to the other profiles at T_{p} = 0.5 s.
Conclusions
This article forms a summary of an extensive set of analyses performed to assess the effect of soil profiles with varying stiffness with depth on the topographic amplification variation across canyon topographies. Topographic amplification factors have been compared between a uniform soil profile and those of varying stiffness considering wavelet input motion imposed in the horizontal direction. It was shown that the position of maximum horizontal topographic amplification A_{hmax} is located between the first and the second natural periods of the 1D soil column transfer function, in agreement with previous studies of ^{Tripe et al. (2013)} and ^{Skiada et al. (2017)}. However, nonuniform soil stiffness profiles mainly increase the topographic amplification factor amplitudes at the canyon crest location in the smaller input motion period range (T_{p} ≤ 0.33 s). Significant differences are observed in the parasitic vertical motion where amplification values could be as high as the horizontal response for some of the nonuniform cases (i.e. Stepped V_{s} profile examined here).