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 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-α 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 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)).

a(t)=βe−αttγsin(2πtTp) (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.

Figure 2 Input motion for T_p= 1 s (α = 2.0, β = 1.5 and γ = 5.0) 

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-held 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.

Ah=PGAtopohPGAffhAh=PGAtopovPGAffh (2)

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., 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 7-02 (2005)}. Both definitions are expressed in equation (3)).

V¯s1=1∑i=1nzi∑i=1nVs,iziV¯s2=∑i=1nzi∑i=1nziVs,i (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.

Table 1 Soil parameters for analyses 

Parameter	Uniform	Linear constant	Toro Vs	Stepped Vs
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
V¯s,m/s	500	493	529	500	250 - z < 62.5 m
					750 - z > 62.5 m
Vs125, 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

Figure 3 Shear and bulk moduli and shear wave velocity profiles with depth 

The soil mass density is equal to ρ = 2.0 Mg/m³ and the horizontal coefficient of earth pressure K₀ = 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₂(ω)| 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₂(ω)| 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)}.

Figure 4 Transfer Function |F₂(ω)| – 1D ground response for the profiles studied 

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₂(ω)| 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₂(ω)| reaches its minimum value, due to the normalisation process by the 1D soil column (^{Skiada et al., 2017}).

Figure 5 Topographic amplification and transfer function for Uniform and Linear Constant profiles 

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.

Figure 6 Maximum normalised response versus predominant period T_p: a) horizontal amplification and b) vertical amplification 

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 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).