Canyon topography effects on ground motion: Assessment of different soil stiffness profiles

David Solans, Evangelia Skiada, Stavroula Kontoe and David M. Potts Department of Civil and Environmental Engineering, Imperial College London, South Kensington Campus, London SW7 2AZ, United Kingdom, david.solans16@imperial.ac.uk, evangelia.skiada11@imperial.ac.uk, stavroula.kontoe@imperial.ac.uk, d.potts@imperial.ac.uk Efectos de topografía de cañón en movimientos sísmicos: Evaluación de diferentes perfiles de rigidez del suelo

The effect of topography on ground motion has been well recognized during numerous earthquakes. Several studies present observational evidence from destructive earthquakes, where the damage is higher in the vicinity of hills and near slope crests. Furthermore, a number of numerical studies aimed to reproduce this phenomenon using different numerical methods, e.g. Finite Elements, Finite Differences and Boundary Elements have been carried out. Most of these investigations involve parametric studies, considering different variables. However, one of the assumptions of these studies is a homogeneous soil stiffness with depth, which is not in most cases realistic. This article investigates the effects of canyon topography on ground motion considering different soil stiffness profiles over a rigid bedrock. Three soil profiles with stiffness variation with depth are examined and compared to the case of a soil layer of uniform stiffness. An additional analysis of a twolayer medium lying above half-space is also considered. Time domain numerical analyses are carried out with the Imperial College Finite Element Program ICFEP, considering linear elastic soil behaviour over rigid bedrock. The input motions are wavelets of harmonic nature, modified by a Saragoni and Hart (1973) temporal filter. These wavelets with a characteristic. pulse period T p in the range of 0.1 s to 2 s are analysed. This study confirms that the topographic amplification extrema are located between the natural periods of the corresponding one-dimensional free-field profile in agreement with recent previous studies. Furthermore, the amplitude of the topographic amplification peaks is shown to change for the different examined soil stiffness profiles. Keywords Saragoni y Hart (1973). Estas ondas con periodos característicos T p en un rango de 0.1 s a 2 s son analizadas. Este estudio confirma que la respuesta normalizada máxima se localiza entre el primer y segundo periodo natural correspondiente a un perfil unidimensional en condición de campo libre de acuerdo a otras investigaciones recientes. Además, la amplitud de la amplificación topográfica muestra modificaciones para los distintos perfiles de suelo examinados.

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) (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 andPapadimitriou, 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. (2004Assimaki et al. ( , 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 half-space 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
Two-dimensional 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, crest-to-crest 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 in-plane shear waves (SV waves) the soil is modelled as a linear visco-elastic material with varying stiffness with depth. Figure 1: Geometry of the domain considered in finite element analysis (Skiada et al., 2018) The study was carried out using the Imperial College Finite Element Program ICFEP (Potts and Zdravković, 1999), employing the generalised-a 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, Dl, was considered as Dl ≤ λ/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 Dl = 5 m was used.
The domain reduction method DRM is used to ensure that free-field 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 two-step 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 a, b and g 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 a = 2.0, b = 1.5 and g = 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 free-field 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 free-field crest acceleration response. The free-field 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 free-field 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 free-field (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 de-coupling 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 free-field ground motion in the case of SV wave incidence. This is because there is no vertical component of the ground motion in the free-field 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 tri-linear 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., 2017Skiada et al., , 2017a. Solans, D., Skiada, E., Kontoe, S. and Potts, D. (2019). Canyon topography effects on ground motion: Assessment of different soil stiffness profiles.
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.  The soil mass density is equal to r = 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, x = 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 (w)| 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 (w)| 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 (w)| 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 (w)| reaches its minimum value, due to the normalisation process by the 1D soil column (Skiada et al., 2017). Solans, D., Skiada, E., Kontoe, S. and Potts, D. (2019). 25, 51-58 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.16-0.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 Figure 6: Maximum normalised response versus predominant period T p : a) horizontal amplification and b) vertical amplification Solans, D., Skiada, E., Kontoe, S. and Potts, D. (2019). Canyon topography effects on ground motion: Assessment of different soil stiffness profiles. 25, 51-58 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, non-uniform 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 non-uniform cases (i.e. Stepped V s profile examined here).