Contents lists available at ScienceDirect

Computers and Geotechnics

journal homepage: www.elsevier.com/locate/compgeo

Research Paper

An energy-based interpretation of sand liquefaction due to vertical ground motion

CrossMark

Vasiliki Tsaparlia'*, Stavroula Kontoea, David M.G. Tabordaa, David M. Potts

a Dep. of Civ. Eng., Imperial College, London, UK

ARTICLE INFO

Article history:

Received 17 December 2016 Received in revised form 2 May 2017 Accepted 3 May 2017

Keywords:

Liquefaction

Vertical ground motion

Energy dissipation

Hysteretic damping ratio

Post-resonance response cycles

ABSTRACT

In several recent earthquakes, high vertical ground accelerations accompanied by liquefaction were observed. Downhole records have also shown that large vertical accelerations do not necessarily originate from the source, but rather get amplified towards the ground surface. Given the advantages of energy-based interpretation of liquefaction triggering due to S-waves, this approach is used together with finite element analyses to investigate vertical motion amplification and ensuing liquefaction. The results show the importance of the post-resonance response cycles, while hysteretic damping based on total stresses, accounting for the water in the pores, is shown to be very low, explaining the observed amplification. © 2017 The Authors. Published by Elsevier Ltd. This is an open access article under the CC BY license (http://

creativecommons.org/licenses/by/4.0/).

1. Introduction

It is widely acknowledged that energy-based concepts are appropriate for evaluating and predicting seismic damage and liquefaction potential: they can implicitly account for the non-linear behaviour of soils, as well as for the distinctive characteristics of the input ground motions, including their irregular and multidirectional nature [1-3].

Davis & Berrill [4] and Berrill & Davis [5] were amongst the first ones to propose an energy-based method (EBM) for liquefaction evaluation on the basis of the assumption introduced by Nemat-Nasser & Shokooh [6], according to which the relationship between unit volume system energy dissipation and the excess pore pressures that develop from the onset of undrained loading and up to liquefaction is unique. The latter assumption has been experimentally proven for both irregular and harmonic loading patterns by a number of researchers over the years [3,7-13] and, as such, more EBMs for liquefaction triggering have been developed [10,11,1321]. Along those lines Green & Terri [1] developed an alternative energy approach to the Palmgren-Miner theory [22,23], originally adopted by Seed et al. [24] for the conversion of an irregular time-history to an equivalent number of uniform stress amplitude

* Corresponding author at: Skempton Building, Imperial College, London SW7 2AZ, UK.

E-mail addresses: vasiliki.tsaparli10@imperial.ac.uk (V. Tsaparli), stavroula. kontoe@imperial.ac.uk (S. Kontoe), d.taborda@imperial.ac.uk (D.M.G. Taborda), d.potts@imperial.ac.uk (D.M. Potts).

cycles. Magnitude scaling factors based on energy concepts for liquefaction evaluation were also proposed by Arango [25] and Green [20], while a number of energy-based pore pressure generation models for sands and silt-sand mixtures have been developed e.g. [4,5,14,26].

More recently, Taborda et al. [27] developed a numerical procedure in finite element analysis for the assessment of dissipated energy and damping ratio in time and space as a result of hysteresis. The implementation of the proposed algorithm is performed in generalised three-dimensional stress space and has uncoupled deviatoric and isotropic mechanisms such that the source that dissipates energy can be isolated and understood.

Given the advantages of energy methods in liquefaction evaluation during S-wave propagation, this paper aims to extend the conclusions on vertical motion amplification and sand liquefaction drawn in an earlier theoretical study by the authors [28], which found resonance to be of fundamental importance, by adopting an energy approach in finite element analyses of P-wave propagation modelling. The motivation behind this emanates from the significantly high vertical ground accelerations recorded during the 22nd February 2011 seismic event in Christchurch, New Zealand, which were accompanied by extensive liquefaction [29,30]. Tsaparli et al. [28] analysed the Fourier spectra of strong motion stations in central Christchurch with the highest registered vertical accelerations and concluded that based on their stratigraphy and ground water conditions, the scenario of resonance in the vertical direction was very likely. As such, the first part of the study presented in this paper focusses on the number of cycles that

http://dx.doi.org/10.1016/j.compgeo.2017.05.006 0266-352X/© 2017 The Authors. Published by Elsevier Ltd.

This is an open access article under the CC BY license (http://creativecommons.org/licenses/by/4.0/).

contribute to soil liquefaction during P-wave propagation using the concept of energy dissipation. The second part provides an explanation for the amplification of vertical seismic motion towards the ground surface, as shown by downhole array field records [31,32], by examining the hysteretic damping ratio induced by both total and effective stress-strain loops. Furthermore, it is noted that this is the first time that the predictions of an advanced bounding surface plasticity model for sand liquefaction, which can accurately replicate a range of features of soil response when subjected to cyclic loading, are interpreted on an energy basis.

2. Vertical ground motion and soil resistance to liquefaction

The possible mechanism of liquefaction triggering due to vertical propagation of compressional waves (P-waves) has been demonstrated using real earthquake records in Tsaparli et al. [28]. In particular, it was shown that resonance in compressional mode for the case of P-wave propagation through a fully saturated sand deposit can lead to significant deviatoric stresses, plasticity and possibly liquefaction. The mechanism is briefly demonstrated herein utilising harmonic waves. For this purpose, a hypothetical fully-saturated uniform 20 m deep level-ground soil deposit with the water table (WT) at ground level, consisting of Fraser River Sand (FRS) of 40% relative density is modelled. The properties of FRS are presented in Tablel.

The maximum shear modulus, Gmax, is given by Eq. (1) [33]:

G b ■ pref ( p' V/2

Gmax -(0.3 + 0.77 ■ e2) UeJ (1)

where B is the elastic shear modulus constant and p¡.ef a reference pressure. Assuming a constant Poisson's ratio, the average compres-sional wave velocity in the deposit is governed by the stiffness of the pores' water and is equal to 1604 m/s, resulting in a fundamental frequency of 20.05 Hz.

Two harmonic waves with a total duration of 10 s are subsequently used to demonstrate resonance and sand liquefaction due to the vertical seismic component. Both have an acceleration amplitude of 0.6 m/s2, however, their predominant frequency differs: the first one has a predominant frequency of 30 Hz (denoted as f30_A0.6) resulting in approximately 300 cycles, whereas in the second one it matches the fundamental frequency of the deposit, i.e. 20 Hz (denoted as f20_A0.6), resulting in 200 cycles. Both acceleration time-histories have been baseline corrected using the computer software SeismoSignal v5.1.0 [34].

All finite element (FE) analyses are non-linear, elasto-plastic, effective stress-based and were carried out with the Imperial College Finite Element Program (ICFEP, [35]), using the u-p hydro-mechanically coupled formulation [36]. This is a simplified form of the complete u-p-w coupled consolidation FE formulation, in which, apart from the solid displacement, u, and the pore water pressure, p, the fluid velocity relative to the solid, w, is an additional degree of freedom. The u-p formulation, however, was identified by Zienkiewicz et al. [36] to be a valid alternative for problems involving frequencies found in common earthquake

Table 1

Material properties for Fraser River Sand.

Properties Value

Gs: specific gravity of soil particles 2.720

e0: initial void ratio 0.812

ysat: saturated bulk unit weight 19.120 kN/m3

K0: earth pressure co-efficient at rest 0.440

v: Poisson's ratio 0.200

k: permeability 4.200E—04 m/s

engineering applications, as in these cases the effect of fluid acceleration relative to the solid as well as its convective terms were not found to be significant. It should be noted that, given the soil permeability and the characteristics of the ground motions and deposits used in this study, the considered problem lies within the range over which the u-p formulation is valid [36]. The FE mesh is a soil column consisting of 8-noded isoparametric quadrilateral elements with pore water pressure degrees of freedom assigned to the four corner nodes (hybrid elements). Plane strain conditions are assumed in all analyses. The size of the elements in the direction of wave propagation was chosen such that it is in agreement with the recommendations by Bathe [37] for the frequency range of the input motions under consideration. Additionally, it should be noted that, when modelling vertical propagation of compres-sional waves in a fully saturated sand deposit, the P-wave velocity is governed by the bulk stiffness of the water, which remains fairly constant. Despite this, in the present work, a conservative approach was adopted and it was assumed that the soil stiffness could be reduced to 20% of its maximum small-strain value. This resulted in element dimensions of 0.1 x 0.1 m2. In terms of boundary conditions, horizontal displacements have been restricted at the bottom boundary of the mesh, while one-dimensional soil response is ensured by tying the displacements in both the horizontal and vertical direction along nodes of the same height [38]. The pore pressures at the nodes along the lateral boundaries are also tied, ensuring one-dimensional fluid flow. The hydraulic regime is additionally defined by a condition of no flow at the base of the mesh and zero pore water pressure change at the top boundary, thus allowing drainage to take place one-dimensionally through the ground surface [35]. The final boundary condition is that of acceleration time-history which is applied incrementally in the vertical direction at the base of the mesh. An accelerated modified Newton-Raphson scheme with a sub-stepping stress point algorithm was employed for the solution of the non-linear constitutive equations [35], while the time-integration scheme is the generalised a-method of Chung & Hulbert [39] with a spectral radius at infinity, pM, of 0.818 [40,41]. In terms of time-step, a value as small as 0.003 s was deemed essential to ensure accuracy for the deposit and input motions under consideration.

To model sand behaviour and liquefaction, a two-surface bounding surface plasticity model has been implemented in ICFEP in three-dimensional stress space [42,43]. This is a modified version of the Papadimitriou & Bouckovalas [44] bounding surface plasticity model, based on the original proposal by Manzari and Dafalias [45]. Details of the modifications as implemented in ICFEP, including a power law for the Critical State Line, an altered expression for the hardening modulus and the introduction of a secondary yield surface, can be found in Taborda [42] and Taborda et al. [43].

The model parameters for Fraser River Sand were established by Klokidi [46] and have been included in Table 2. The meaning of these is explained in detail in Taborda [42] and Taborda et al. [43] and is not repeated herein for brevity. The performance of the model at element level in replicating the response of FRS is included in Appendix A. The stiffness degradation and damping variation curves at a mean effective stress level of 60 kPa, corresponding approximately to the mid-depth of the 20 m deep FRS deposit, have also been reproduced and included in Appendix A. The computed curves are compared with empirical curves from the literature [47]. It should be noted that the FRS was chosen because of an abundance of element testing readily available (see Appendix A).

Fig. 1 shows the input and ground surface acceleration time-histories and Fourier Spectra for each of the two motions considered, while Fig. 2 depicts the corresponding initial and final mean

Table 2

Model parameters for Fraser River Sand [46].

Model Value Model Value Model Value

Pi-ef (kPa) 100.00 A0 1.00 he 0.119

(eCs)ref 0.95 m 0.065 c 1.016

k 0.05 PYS (kPa) 1.00 emax 0.97

n 0.60 B 422.00 a 1.00

mc 1.38 ai 0.44 b 0.00

Me 1.00 j 2.00 l 1.00

kb 2.50 Yi 7.95E—04 Ha 14000.00

kd 1.80 V 0.20 f 1.16

•Critical state line: eCS = (eCS)ref -

Fig. 1. Ground surface and input acceleration time-history of (a) f30_A0.6, (c) f20_A0.6 and ground surface and input Fourier spectra of (b) f30_A0.6 and (d) f20_A0.6.

0 40 p'(kPa) 80 120 0 40 p'(kPa) 80 120

Fig. 2. Computed initial and final mean effective stress profiles during the strong motion (a) harmonic motion with a predominant frequency of 30 Hz (f30_A0.6) (b) harmonic motion with a predominant frequency of 20 Hz (f20_A0.6).

effective stress profiles along the 20 m deep deposit. Clearly, for f30_A0.6, as the amplitudes of the components of the input motion close to the fundamental frequency of the deposit (i.e. 20 Hz) are zero, the resulting ground surface motion is very similar to the input one (Fig. 1a). As a result, the response is elastic with no change in mean effective stresses (Fig. 2a). Conversely, in the case of f20_A0.6, as the predominant frequency coincides with the natural frequency of the system, resonance leads to significant amplification towards the ground surface (Fig. 1c and d), leading to the development of substantial deviatoric stresses and plastic strains. It is evident from Fig. 2b that the seismic loading is now suffi-

ciently strong to trigger liquefaction, even though the deposit has been subjected only to vertical ground motion.

3. Number of cycles in liquefaction triggering due to vertical seismic motion

The above findings indicate that the peak ground acceleration of the input motion is not of significant importance when it comes to vertical motion and cyclic strength. Similarly, it would seem that for the phenomena described herein the number of cycles of the input motion is not fundamental: input motion f30_A0.6 is charac-

terised by a larger number of cycles compared to f20_A0.6, however, it is the latter that leads to soil plasticity and liquefaction. As such, to better understand the fundamental mechanism of liquefaction due to the vertical seismic component, two distinct aspects of the number of cycles are investigated further in the present study: the number of post-resonance response cycles and the number of cycles of the input motion.

3.1. Number of post-resonance response cycles

In Tsaparli et al. [28] two different earthquake records were used to highlight the importance of the frequency content of the input vertical motion: the outcrop 22nd February 2011 Christch-urch motion [48] as recorded at the Lyttelton Port Company (LPCC) strong motion station (SMS) and the 20th May 1986 Lotung, Taiwan, seismic motion as recorded at 47 m depth in a downhole array [49]. The modelled duration of the two motions were 20 s and 40 s, respectively, while the Lotung vertical component was scaled up to have the same peak ground acceleration (PGA) with that of the Christchurch motion (i.e. 0.4 g). It was then shown that the Christchurch vertical component (CV) is characterised by a much wider frequency content than the Lotung one (LV), with frequencies up to 25 Hz. As a result, when CV was used as input excitation in a 40 m deep FRS deposit (CV_40), with a fundamental frequency for P-waves equal to 10.3 Hz, resonance led to liquefaction across the entire depth after 20 s of strong motion. Conversely, when the LV motion was modelled in the 40 m deep deposit (LV_40), despite the double duration of strong motion compared to CV_40 (i.e. 40 s versus 20 s), the response was elastic as resonance did not take place. For this reason, a deeper deposit of 166 m was used to model the propagation of LV (LV_166), in order to shift the natural frequency to 2.46 Hz and hence to a range of frequencies predominant in the LV input motion. The result was the triggering of liquefaction due to resonance now taking place. Cyclic stress ratio time-histories were also calculated for the various analyses, with the computed maximum amplitudes justifying the observed response. These analyses and conclusions have been exhaustively discussed in the relevant study and the reader can refer to that for more details [28].

In the above analyses the assumption of a constant relative density of 40% results in a variable state parameter, W [50], and, thus, a variation in the response of the sand with depth across the deposit: the state parameter ranges from a value of -0.138 at 0 m depth to -0.055 at 40 m and 0.000 at 93 m depth, increasing to a value of 0.057 at the base of the 166 m deep FRS deposit.

It should be noted that the minimum modelled values of the state parameter correspond to less denser-than-critical states (i.e. the state parameter is less negative) than those estimated at similar depths at a number of strong motion stations in Christchurch, New Zealand, which exhibited surface manifestation of liquefaction, as well as cyclic mobility characteristics in the recorded ground surface acceleration time-histories. A characteristic example is the Pages Road Pumping (PRPC) station which, apart from showing evidence of liquefaction [51], also recorded the second largest peak vertical acceleration, equal to 1.63g, during the Mw 6.2 February Christchurch seismic event [48]. The field state parameter at the PRPC station, as established through a CPT-W correlation by Robertson [52], attained values between -0.283 and -0.097 at depths ranging from 3 to 20 m, where the sand layer lies, with W values more negative than -0.2 (i.e. corresponding to clear denser-than-critical states) across the first 7 m of the sand stratum. Boreholes and CPT's in the vicinity of the PRPC station were obtained through the New Zealand Geotechnical Database [53]. It is evident that the range of the state parameter in the deposit modelled herein corresponds to a less denser-than-critical state than

that at the PRPC station between 3 and 20 m depth (a simulated W value of -0.120 at 3 m depth to a value of -0.083 at 20 m depth).

The conclusions made above are further extended herein by looking at the number of post-resonance response cycles, which essentially represent the duration of motion once resonance has taken place. Motivation for this investigation was given by the fact that, although liquefaction took place in analysis LV_166 towards the end of the 40 s duration, when the Christchurch motion duration is used as the basis for comparison (i.e. 20 s), plasticity is triggered in LV_166, but not liquefaction. This is clearly shown in Fig. 3a where the mean effective stress profiles in the 166 m deep deposit have been plotted for the Lotung motion at 20 s and 40 s, along with the initial profile.

To understand the key mechanism behind this, the acceleration time-history at ground surface for analysis LV_166 has been plotted in Fig. 3b and is compared with the respective acceleration time-history from the CV_40 analysis which, in the duration of 20 s, led to liquefaction. A key distinction between the two time-histories in the first 20 s is the smaller number of response cycles in the case of analysis LV_166. This is not surprising if one considers that a deeper deposit is characterised by a larger natural period and hence would lead to a smaller number of response cycles compared to a shorter depth. This could present a plausible explanation for the non-occurrence of liquefaction in the first 20 s in LV_166 and, similar to S-wave propagation, shows the prominent role of the response cycles in vertical motion and liquefaction triggering.

Additionally, it is interesting to note that in the case of LV_166 resonance seems to take place only after about 9 s of strong motion. To further examine this, the Fourier spectrum of LV from 0 to 9 s and from 9 to 19 s has been separately plotted in Fig. 4. The fundamental frequency (f1) of the 166 m deep FRS deposit for P-waves at 2.46 Hz is also shown for clarity. It is clear that the first 9 s of motion consist of practically no components close to f1, with most of them incorporated within the 9-19 s duration.

To verify whether liquefaction would be triggered within 20 s in LV_166 if resonance would have taken place from the beginning of the strong motion (i.e. similar to CV_40), the Lotung vertical component was modified by removing the first 0-9 s and combing the 9-19 s part twice in a total duration of 20 s (herein designated as LV_II). The resulting acceleration time-history and corresponding Fourier spectrum are shown in Fig. 5. This artificial motion was also baseline corrected [34] prior to its use in the finite element analyses.

As anticipated, Fig. 6a shows that resonance in the case of LV_II_166 takes place from the beginning of the strong motion, resulting in liquefaction towards the end of the 20 s duration (Fig. 6b). The effective stress paths in p'-J space, where J is the second deviatoric stress invariant, at mid-depth of the deposit during the first 20 s of motion for the two analyses under consideration have also been included in Fig. 7. In both cases, as resonance takes place, substantial deviatoric stresses develop due to the changes in the direct effective stresses, leading to a reduction in the mean effective stress, with the stress path moving towards the isotropic axis and the origin of space. As expected, the paths are initially identical. However, as resonance takes place only in the second half of the 20 s duration in LV_166, the number of post-resonance response cycles is not sufficient to lead to liquefaction and zero effective stresses. This proves not only the importance of resonance, as previously identified [28], but also that of the number of post-resonance response cycles in liquefaction triggering due to the vertical motion. This establishes a fundamental contrast with problems involving horizontal motions since, in such cases, resonance may aid the triggering of liquefaction, but it is not a prerequisite.

The conclusions made above can be better understood by studying the energy dissipated in the above analyses. The dissipated

Fig. 3. (a) Mean effective stress profile at 0 s, 20 s and 40 s of strong motion and (b) ground surface acceleration time-history as registered for the Lotung seismic event in the 166 m deep FRS deposit.

Fig. 5. Acceleration time-history and respective Fourier spectrum of the modified input Lotung vertical component (LV_II).

0 400 P'(kPa) 800

0 5 10 15 20

t (s) (b) (a)

Fig. 6. (a) Comparison of the computed ground surface acceleration time-history and (b) mean effective stress profiles at 20 s for the vertical component of the modified Lotung seismic event (LV_II_166) and the original analysis (LV_166).

energy time-histories have been based on vertical effective stress- Taborda et al. [27]. The simplification for one-dimensional loading vertical strain loops and have been calculated using the algorithm conditions of the general methodology for the calculation of con-

implemented in ICFEP in three-dimensional stress space by tinuous histories of energy dissipation and damping is shown

100 ■

:- LV 166 20 s LV II 166 20 s

83 m depth -i—i—1—i—i—i—i—i—

p' (kPa)

Fig. 7. Computed stress paths in mean effective - deviatoric stress space for the first 20 s of ground motion at mid-depth of the 166 m deep FRS deposit for the modified Lotung seismic event (LV_II_166) and the original analysis (LV_166).

schematically in Fig. 8 and is applicable to irregular motions where the stress-strain loops are not necessarily closed.

The accumulated energy, Eacci, between the last reversal point and the current increment i in the analysis is represented by the grey area enclosed by the stress-strain curve, and is given by Eq. (2):

Eacc,i —

rv i rv rev ) ' dev

where £v>rev and r, rev are the vertical strain and vertical effective stress corresponding to the last reversal point, respectively, while £v>i and r, i are the corresponding current values. A reversal point is registered whenever a strain reduction takes place. Similarly, the elastic energy of the material between these two points is represented by the hatched area in Fig. 8 and is given by Eq. (3):

Eel i — 2'

(rV,i rV,rev) ' (ev,i ev,rev)

In this way, if it is assumed that increment i is the subsequent reversal point, then the dissipated energy and elastic energy of the hypothetical full loop between the last reversal and increment i are given by Eqs. (4) and (5):

Fig. 8. Definition of accumulated and elastic energy based on a stress-strain curve (adapted from Taborda et al. [27]).

Ediss , — Aloop — 2 ' (Eacc,i Eel ,i) Eel i

Eel — Ael —

Having defined the above quantities, a continuous evaluation of the damping ratio can also be performed by combining Eqs. (2) and (3) in the following form:

1 Ediss _ 2 ' (Eacc,i - Eel,i)

4 ' p Eel p ' Eel i

For more details on the above algorithm and its implementation for general stress space, the reader can refer to Taborda et al. [27].

The dissipated energy time-histories for analyses LV_166 and LV_II_166 are shown in Fig. 9 for two depths in the deposit: 30 m and 83 m. It is clear that practically no energy is dissipated in LV_166 during the first 10 s of strong motion, while a continuous dissipation of energy takes place from the beginning of the motion in LV_II_166, justifying the larger amounts of plasticity and liquefaction that occur in the latter. A similar rate of energy dissipation to that between 9 and 19 s in LV_166 is noted for the first 10 s in LV_II_166, although somewhat lower, as the amplitudes of acceleration in the latter analysis are smaller when compared to those of the former once resonance has taken place. This has been shown in Fig. 6a and can be explained by the slightly smaller input Fourier amplitudes close to f1 in the modified Lotung vertical component in the first 10 s compared to the original motion for the full duration. The latter can be seen in Fig. 4. Additionally, the rate of energy dissipation from 10 to 20 s in LV_II_166 is substantially different from that during 0 to 10 s, reflecting the effect of the stressstrain state on the energy dissipation: at 10 s the stress-strain state is significantly altered from that at 0 s as permanent excess pore water pressures and plastic strains have developed.

Finally, it is of interest to verify the consistency of the above observations with the computed mean effective stress-ratio variation with various quantities (i.e. time, dissipated energy, number of response cycles, representative cycles based on stress and representative cycles based on energy) for the above two analyses. The mean effective stress ratio, rp, has been defined in Tsaparli et al. [28] according to Eq. (7):

r _ ap'

rp — —7~

where Dp' is the change in mean effective stress since the start of the analysis and p'0 is the mean effective stress level after consolidation, prior to the application of the dynamic loading. rp replaces the commonly used excess pore water pressure ratio, ru, defined as the ratio of the change in pore water pressure, Du, over the vertical effective stress after consolidation, r,0, as the former is unaffected by total stress changes which can be substantial during P-wave propagation. These total stress variations are transmitted to the water phase and can lead to ru values greater than unity even when liquefaction has not taken place. Liquefaction occurrence is defined here by an rp value equal to about 0.9, with a value of unity corresponding to complete loss of soil strength (initial liquefaction [54]), similar to ru.

Fig. 10 shows rp variations against time (t), dissipated energy (Ediss), number of response cycles (Ntot), representative cycles based on stress (Nstressrepr) and representative cycles based on energy (Nenergy repr), as registered during the 20 s of strong motion for analyses LV_166 and LV_II_166. Two depths in the 166 m deep FRS deposit have been selected: 30 m and 83 m. From the above quantities, the dissipated energy has already been defined according to Eq. (4). A loading cycle or loop and its corresponding dissipated and elastic energies, Ediss,loop and Eel,loop, are calculated and stored as described earlier (Eqs. (4) and (5)) and schematically

ev,rev

shown in Fig. 8, whenever two consecutive reversals, rev; and revi+i, are registered. As such, the hysteretic damping ratio of a cycle is given by Eq. (8):

■diss,loop

4- p E,

el,loop

The time-history of the response cycles is then obtained by assuming that the registration of two consecutive reversals corresponds to half a cycle.

To calculate the number of representative cycles based on stress or energy, the stress amplitude of each loop is also calculated as half the stress range between two consecutive reversals. Only the loops which are characterised by a stress amplitude or dissipated energy, Ediss loop, larger than or equal to 30% of the maximum stress amplitude or Edmisasx,loop registered during the whole loading history are considered. These are essentially introduced to account for the post-resonance response cycles. The 30% factor was chosen to be consistent with the Seed et al. [24] methodology, which converts an irregular stress or acceleration time-history to an equivalent number of uniform amplitude cycles. In this approach, cycles of stresses or acceleration with amplitude less than 30% of the maximum value are neglected, as the corresponding weighting factors are insignificant in terms of liquefaction triggering.

As expected, the results in terms of mean effective stress ratio time-histories depicted in Fig. 10 show that rp increases from the early stages of seismic loading in LV_11_166, reaching values of 1 towards the end of the 20 s duration, which indicates the occurrence of liquefaction. Conversely, the value of rp in LV_166 starts increasing only after about 9 s of motion, as no resonance has taken place up to this point. Between 0 s and 9 s no permanent changes in stresses develop with only high frequency oscillations being present in the rp time-histories, due to the total stress changes and the assumption of compressible pore fluid. The maximum value that rp attains in this case is 0.6 and 0.4 for the depths of 30 m and 83 m, respectively, indicating plasticity but no liquefaction.

When rp is plotted against dissipated energy, E^iss, it becomes clear that the results from both analyses initially coincide, however, as anticipated, more energy is dissipated overall in the case of LV_11_166. The representation of rp against the response cycles, Ntot, also proves that the amplitude of cycles during the first 9 s of strong motion in analysis LV_166 is not significant as no permanent change in mean effective stresses, p', takes place. Finally, when rp is plotted against the number of representative cycles based on stress, Nstress,repr, or energy, Nenergy,repr, the variations from the two analyses once again initially are very similar. Analysis LV_11_166, however, is shown to result in almost double the number of representative cycles when compared to LV_166. This indicates that the selection of a factor of 30% in the definition of

representative cycles is adequate for the problem under consideration, especially for the case of those based on energy.

3.2. Number of cycles of the input motion

1t emanates from the conclusions drawn above that, in the case of vertical seismic loading and liquefaction, the number of cycles of the input ground motion are not of fundamental importance. To demonstrate this for the case of real earthquake records, the vertical component of the 20th September 1999 Chi-Chi, Taiwan, seismic event with a moment magnitude of Mw = 7.7 [55], as recorded at 52 m depth in a downhole array at the Hualien site, was used to carry out an additional site response finite element analysis in 1CFEP, utilising the 40 m deep FRS deposit. The aim of this analysis (denoted as ChiV_40) is to compare it with the results of the original CV_40 analysis for the Christchurch motion, in which plasticity and liquefaction took place.

Fig. 11 shows the baseline corrected acceleration time-history and corresponding Fourier spectrum of the Chi-Chi vertical motion (Chi_V). The peak vertical acceleration of the original motion is equal to 0.02g, however, the motion was scaled up to reach a peak acceleration equal to that of the Christchurch vertical component and thus, isolate the effects of frequency content and number of cycles of the input motion. 1t is clear in Fig. 11 that, unlike the Christchurch motion, the chosen Chi-Chi motion does not possess any significant components close to the fundamental frequency of the 40 m deep deposit for P-waves, i.e. 10.3 Hz.

Additionally, the equivalent number of uniform amplitude cycles of the input vertical component of the Christchurch and the Chi-Chi ground motions was calculated in accordance with the Seed et al. [24] methodology. From this it was found that the Christchurch vertical component is equivalent to 4.5 cycles of uniform amplitude at 65% of the maximum acceleration, amax, whereas the Chi-Chi one is equivalent to 18.5 cycles with an amplitude of 65% of amax, reflecting the larger duration of the motion compared to CV.

The similarity of the frequency content of the Chi-Chi motion to that of the Lotung vertical component (i.e. no significant input components at the fundamental frequency of the 40 m deep FRS deposit for compressional waves), as well as the larger number of equivalent uniform amplitude cycles compared to those of the Christchurch vertical component renders Chi_V ideal for studying the effect of the number of cycles of the input motion on liquefaction due to vertical motion.

The results of analysis ChiV_40 are shown in Fig. 12 in terms of input and ground surface acceleration time-histories, as well as initial and final mean effective stress profiles. As there is no resonance for the ChiV_40 case, the ground surface acceleration time-history is similar to the input one and as a result the developed deviatoric stresses are not large enough to lead to plastic strains.

Fig. 10. rp variation with time (t), dissipated energy (Ediss), number of response cycles (Ntot), representative number of cycles based on stress (Nstress, repr) & representative number of cycles based on energy (Nenergy, repr) for the modified Lotung vertical component (LVJL166) and the original one (LV_166) at 30 m & 83 m depth.

Fig. 12. (a) Input and computed ground surface acceleration time-histories for the vertical component of the Chi-Chi seismic event and (b) initial and final mean effective stress profiles during the strong motion as computed for the vertical component of the Chi-Chi seismic event in the 40 m deep FRS deposit (ChiV_40).

Fig. 13. Dissipated energy time-histories as computed for the Chi-Chi (ChiV_40) and Christchurch (CV_40) vertical component analyses at 10 m and 30 m depth in the 40 m deep FRS deposit.

Consequently, the mean effective stress profile remains unaltered. The results in terms of dissipated energy time-history at 10 m and 30 m depth, shown in Fig. 13, prove that no significant hysteresis takes place as the response is elastic and very small amounts of energy are dissipated.

the transfer function (TF) of CV_40, defined as the ratio of the ground surface Fourier Spectrum over the input one, with the analytical TF for the steady state solution of a harmonic wave propagating through visco-elastic soil on rigid rock [56]. As it can be

4. Soil amplification of the vertical ground motion

The final aspect investigated, in terms of energy and damping, is the mechanism of amplification of vertical motion towards the ground surface. Several field records of downhole arrays [31,32] have shown that, even when liquefaction takes place, the vertical motion is amplified towards the ground surface.

Tsaparli et al. [28] showed that the ground surface acceleration of the Christchurch vertical component in analysis CV_40 retained high amplitudes even after the first 10 s, when the strong part of the input motion is completed. This can be seen here in Fig. 3b. 1t was initially postulated that this was a result of low levels of damping, resulting from the non-hysteretic behaviour of the water in the pores which governs the response in the vertical direction [28]. This damping can be estimated by matching the peaks of

0 2 4 6 8 10 12 14

f (Hz)

Fig. 14. Comparison of the predicted transfer function for the Christchurch vertical component analysis of the 40 m deep FRS deposit (CV_40) and the analytical transfer function for the steady-state solution of the propagation of a harmonic wave through visco-elastic soil overlying rigid bedrock.

seen in Fig. 14, a value of damping as small as 0.4% is required to provide a good match.

In the current study, the mechanism of amplification of the vertical motion, which takes place even in the case of significant plasticity and occurrence of liquefaction, is further explored by computing the hysteretic deviatoric damping ratio time-histories at various depths within the 40 m deep FRS deposit for the analysis using the Christchurch vertical component (CV). These have been calculated based on Eq. (6) for continuous damping ratio time-histories described earlier [27]. Two different values of damping ratio have been calculated:

• The first one, denoted as D1, is the effective stress mechanism based on vertical effective stress-vertical strain (o^ - ev) loops; this reflects the solid skeleton.

• The second one, denoted as D2, is the total stress mechanism based on vertical total stress-vertical strain (rv - ev) loops; this accounts for the presence of the non-hysteretic water in the soil's pores.

Fig. 15 compares the above damping ratio time-histories at two locations within the 40 m deep FRS deposit: at 10 m and 30 m depth. As expected, a rapid variation in this quantity is evident, due to the large number of cycles and the fact that the calculated damping ratio is a function of the stress-strain curve between the current point and the last reversal (see Eq. (6)). It can be seen that the damping ratio based on effective stresses (D1) increases to significantly large values, which is consistent with the plasticity registered in the analysis. Tsaparli et al. [28] showed that, at most depths, initial liquefaction in CV_40 took place after about 10 s, towards the end of the strong motion duration. This explains the substantial development of D1 after about 10 s shown in Fig. 15. Conversely, the computed damping ratio based on total stress components (D2) retains significantly low values justifying the amplification seen towards the ground surface. Average values of D2 were found to be equal to 0.47% and 0.27% at 10 m and 30 m depth, respectively, showing a very good agreement with the damping ratio calculated from the analytical TF (Fig. 14). An increasing trend of D2 is obtained for both depths after 10 s of strong motion, justifying the small drop in the acceleration amplitudes from 10 s onwards, as seen in Fig. 3b.

The above observations imply that, if the design earthquake, the stratigraphy details and the fundamental frequency of the deposit under consideration are known, a simple linear elastic frequency domain type analysis with practically zero or very small (<0.5%) viscous (i.e. Rayleigh) damping could predict the potential development of large vertical inertial forces which could lead to com-pressional structural damage, as this has been observed in a number of field observations [57-60]. The compatibility and appropriateness of the conventional frequency domain site response analysis for single-phase P-wave propagation modelling

has been demonstrated by Han [61]. It should be noted that care should be exercised in cases of permeability values higher than 1.0 x 10-3 m/s, as in these cases viscous damping resulting from the interaction between the solid and fluid phases can be quite substantial, reaching values of up to about 10% [61]. This can lead to substantial de-amplification of the vertical acceleration amplitudes, meaning that such small values of damping (<1%) are not suitable for these cases.

5. Conclusions

The present study extends the findings of Tsaparli et al. [28] regarding the occurrence of liquefaction solely due to vertical ground motion propagation at resonance conditions by applying energy-related considerations to the results of a series of finite element analyses. In particular, the role of the number of post-resonance response cycles and of the number of cycles of the input ground motion in triggering liquefaction during P-wave propagation have been analysed by studying time variations of dissipated energy and mean effective stress ratio, as well as by analysing the variation of mean effective stress with number of cycles at various depths across the modelled deposit.

Additionally, the amplification of the vertical motion towards the ground surface, even when liquefaction takes place, is exhaustively discussed by studying two different hysteretic damping ratio time-histories: one based on effective vertical stress-vertical strain loops characterising the solid skeleton and one based on the corresponding total quantities reflecting the presence of water in the soil's pores. From these the following conclusions can be inferred:

In the case where liquefaction is induced by the vertical seismic component, the number of post-resonance response cycles is a key factor, as significant levels of energy dissipate only once resonance takes place. The occurrence of liquefaction in this case will be governed by the number of post-resonance response cycles and their amplitude, as if not sufficient, plasticity can be invoked but no liquefaction. This is significant as the number of response cycles rather than resonance is a prerequisite in the case of S-wave propagation and triggering of liquefaction in a level-ground fully saturated sand deposit.

In a similar manner, the number of cycles in the input motion is shown to be of no significant importance.

When liquefaction takes place due to P-wave propagation, damping as a result of effective stresses increases substantially, in agreement with the plasticity observed. Conversely, the damping calculated based on the variation of total stresses retains very low values, as a result of the non-hysteretic behaviour of the water in the pores, justifying the amplification of vertical motion towards the ground surface even at liquefied sites.

The above implies that the potential of large vertical ground acceleration development and hence compressional structural

damage can be predicted by carrying out a simple linear elastic frequency domain type analysis using appropriate levels of viscous (i.e. Rayleigh) damping, provided that the design earthquake, the stratigraphy details and the fundamental frequency of the deposit are known.

Acknowledgments

The first author would like to gratefully acknowledge the financial support by the Engineering and Physical Sciences Research Council (EPSRC, 1386368). The contribution from the 1nstitute of Earth Sciences, Academia Sinica, is also acknowledged for providing the ChiChi ground motion.

Appendix A

The bounding surface plasticity model has been calibrated for FRS by Klokidi [46] based on monotonic triaxial tests conducted

by [62,63], as well as cyclic undrained direct simple shear tests of [64]. The performance of the model under both monotonic and cyclic solicitations and for different initial states is demonstrated below.

Figs. A.1 and A.2 show a comparison of the computed response with the laboratory data under monotonic triaxial conditions for different mean effective stress levels after consolidation and at distinct relative densities, respectively. The response is shown in terms of stress-strain curves, effective stress paths and excess pore water pressure generation curves. 1t is shown that the model, through the state parameter concept, is able to capture successfully the mechanical behaviour of FRS under different initial states. 1t is also worth noting that the response is accurately simulated for both triaxial compression (Fig. A.1) and triaxial extension (Fig. A.2) solicitations.

The undrained cyclic response of FRS consolidated to a vertical effective stress of 100 and 200 kPa in terms of effective stress path in vertical effective stress-shear stress space and of excess pore

Fig. A.1. Undrained triaxial compression test simulations on FRS with a relative density of 45% isotropically consolidated under different stress levels - (a) deviatoric stressaxial strain curve, (b) stress path in mean effective stress-deviatoric stress space and (c) excess pore water pressure generation against axial strain.

Fig. A.2. Undrained triaxial extension test simulations on FRS of different relative densities isotropically consolidated under a mean effective stress level of 100 kPa - (a) deviatoric stress-axial strain curve, (b) stress path in mean effective stress-deviatoric stress space and (c) excess pore water pressure generation against axial strain.

Fig. A.3. Model performance in undrained cyclic simple shear test simulations on FRS with a relative density of 40% under an initial vertical effective stress of (a) 100 kPa and (b) 200 kPa.

Fig. A.4. Reproduction of drained cyclic triaxial tests on Fraser River Sand - secant shear modulus and damping variation curves versus shear strain amplitude and comparison with empirical curves from the literature.

water pressure development with the number of cycles of load application is illustrated in Fig. A.3 for a relative density of 40%, as the one simulated in the considered boundary value problem. Similar to the monotonic solicitations, the model is shown to perform well under cyclic loading for different initial states.

Finally, the normalised secant stiffness degradation and damping variation curves for FRS corresponding to a mean effective stress level of 60 kPa and a Dr of 40% are shown in Fig. A.4. The curves were determined by interpreting the stress-strain loops obtained from strain-controlled drained cyclic triaxial test simulations. The normalised stiffness and damping ratio values have been plotted against the cyclic shear strain amplitude, yc, where the triaxial shear strain is obtained by subtracting the horizontal strain component from the vertical (axial) strain. The dynamic properties have been estimated based on the behaviour calculated for the 10th stress-strain loop, commonly used for earthquake engineering applications. These are compared against the empirical curves of Darendeli [47]. The latter have been calculated for similar conditions to those characterising the FRS deposit used in the present study (p0 = 60 kPa, OCR = 1, plasticity index = 0%, 10 loading cycles, a loading frequency of 1 Hz). A fairly good agreement can be seen throughout the whole strain range for the normalised stiffness degradation, with the onset of plastic response inferred by the presence of a kink at a strain amplitude of about 5.0E-3%. Reasonable damping values are reproduced by the model in the medium strain range, though the adoption of the basic Masing rules [65]

within the yield surface for the reproduction of hysteresis [44] results in an underestimation of the simulated damping ratio at the very small strain range.

References

[1] Green RA, Terri GA. Number of equivalent cycles concept for liquefaction evaluations—revisited. J Geotech Geoenviron Eng 2005;131:477-88. http://dx. doi.org/10.1061/(ASCE)1090-0241(2005)131:4(477).

[2] Kokusho T. Liquefaction potential evaluations: energy-based method versus stress-based method. Can Geotech J 2013;50:1088-99.

[3] Azeiteiro RJN, Coelho PALF, Taborda DMG, Grazina JCD. Energy-based evaluation of liquefaction potential under non-uniform cyclic loading. Soil Dyn Earthq Eng 2016;92:650-65. http://dx.doi.org/10.1016/ j.soildyn.2016.11.005.

[4] Davis RO, Berrill JB. Energy dissipation and seismic liquefaction in sands. Earthq Eng Struct Dyn 1982;10:59-68.

[5] Berrill JB, Davis RO. Energy dissipation and seismic liquefaction of sands: revised model. Soils Found 1985;25:106-18.

[6] Nemat-Nasser S, Shokooh A. A unified approach to densification and liquefaction of cohesionless sand in cyclic shearing. Can Geotech J 1979;16:659-78.

[7] Simcock KJ, Davis RO, Berrill JB, Mullenger G. Cyclic triaxial tests with continuous measurement of dissipated energy. Geotech Test J 1983;6:35-9.

[8] Towhata I, Ishihara K. Shear work and pore water pressure in undrained shear. Soils Found 1985;25:73-85.

[9] Yanagisawa E, Sugano T. Undrained shear behaviors of sand in view of shear work. In: Proc. 13th Int. Conf. Soil Mech. Found. Eng. (Special Vol. Perform. Gr. Soil Struct. Dur. Earthquakes). New Delhi, India: Balkema Publishers; 1994. p. 155-8.

[10] Figueroa JL, Saada AS, Liang L, Dahisaria NM. Evaluation of soil liquefaction by energy principles. J Geotech Eng ASCE 1994;120:1554-69.

Liang L, Figueroa JL, Saada AS. Liquefaction under random loading: unit energy [39 approach. J Geotech Eng 1995;121:776-81. http://dx.doi.org/10.1061/(ASCE) 0733-9410(1995)121:11(776).

Polito C, Green RA, Dillon E, Sohn C. Effect of load shape on relationship [40 between dissipated energy and residual excess pore pressure generation in cyclic triaxial tests. Can GeotechJ 2013;50:1118-28.

Kokusho T, Mimori Y, Kaneko Y. Energy-based liquefaction potential [41 evaluation and its application to a case history. In: Proc. 6th Int. Conf. Earthq. Geotech. Eng., 1-4 November, Christchurch, New Zealand: 2015. Law KT, Cao YL, He GN. An energy approach for assessing seismic liquefaction [42 potential. Can GeotechJ 1990;27:320-9.

Cao YL, Law KT. Energy approach for liquefaction of sandy and clay silts. In: [43 Proc. 2nd Intern. Conf. Recent Adv. Geotech. Earthq. Eng. Soil Dyn. Pap. No 3.38. Missouri University of Science and Technology; 1991. p. 491-7. Alkahatib M. Liquefaction assessment by strain energy approach. Wayne State [44 University; 1994.

Trifunac MD. Empirical criteria for liquefaction in sands via standard penetration tests and seismic wave energy. Soil Dyn Earthq Eng [45 1995;14:419-26. http://dx.doi.org/10.1016/0267-7261(95)00016-N. Running DL. An energy-based model for soil liquefaction PhD [46 thesis. Washington State University; 1996.

Kayen RE, Mitchell JK. Assessment of liquefaction potential during earthquakes [47

by Arias intensity. J Geotech Geoenviron Eng 1997;123:1162-74.

Green RA. Energy-based evaluation and remediation of liquefiable soils PhD [48

thesis. Virginia Polytechnic Institute and State University; 2001. http://dx.doi. [49

org/10.1061/40744(154)191.

Jafarian Y, Towhata I, Baziar MH, Noorzad A, Bahmanpour A. Strain energy [50 based evaluation of liquefaction and residual pore water pressure in sands using cyclic torsional shear experiments. Soil Dyn Earthq Eng 2012;35:13-28. [51 Palmgren A. Durability of ball bearings-Die lebensdauer von kugella geru. ZVDI 1924;68:339-41.

Miner MA. Cumulative damage in failure. Trans ASME 1945;67:A159-64.

Seed HB, Idriss IM, Makdisi F, Banerjee N. Representation of irregular stress [52

time histories by equivalent uniform stress series in liquefaction analyses. In:

Earthquake Engineering Research Center, College of Engineering, Rep. No.

EERC 75-29. University of California; 1975. [53

Arango I. Magnitude scaling factors for soil liquefaction evaluations. J Geotech

Eng 1996;122:929-36. [54

Green RA, Mitchell JK, Polito CP. An energy-based excess pore pressure

generation model for cohesionless soils. In: John Booker Meml. Symp.. Sydney, [55

Australia: Balkema Publishers; 2000. p. 16-7.

Taborda DMG, Potts DM, Zdravkovic L. On the assessment of energy dissipated

through hysteresis in finite element analysis. Comput Geotech [56

2016;71:180-94. http://dx.doi.org/10.1016/j.compgeo.2015.09.001.

Tsaparli V, Kontoe S, Taborda DMG, Potts DM. Vertical ground motion and its [57

effects on liquefaction resistance of fully saturated sand deposits. Proc R Soc A

Math Phys Eng Sci 2016:472. http://dx.doi.org/10.1098/rspa.2016.0434.

Bradley BA. Recorded ground motions from the 22 February Christchurch [58

earthquake. Taormina, Italy: Second Int. Conf. Performance-Based Des. Earthq.

Geotech. Eng.; 2012.

Lee RL, Franklin MJ, Bradley BA. Characteristics of vertical ground motions in [59 the Canterbury earthquakes. New Zeal. Soc. Earthq. Eng. Annu. Conf.. Wellington: University of Canterbury; 2013.

Yang J, Sato T. Interpretation of seismic vertical amplification observed at an array site. Bull Seismol Soc Am 2000;90:275-85. [60

Yang J, Sato T, Savidis S, Li XS. Horizontal and vertical components of earthquake ground motions at liquefiable sites. Soil Dyn Earthq Eng 2002;22:229-40. [61

Hardin BO. The nature of stress-strain behaviour for soils. Proc. ASCE Geotech. Eng. Div. Spec. Conf., vol. 1, Pasadena, CA.; 1978, p. 3-90.

Seismosoft. Seismosignal - a computer program for signal processing of [62 strong-motion data; 2014.

Potts DM, Zdravkovic L. Finite element analysis in geotechnical engineering: [63 theory. London: Thomas Telford; 1999.

Zienkiewicz OC, Chang CT, Bettess P. Drained, undrained, consolidating and dynamic behaviour assumptions in soils. Géotechnique 1980;30:385-95. [64 http://dx.doi.org/10.1680/geot.1980.30.4.385.

Bathe KJ. Finite element procedures. New Jersey: Prentice Hall; 1996. Zienkiewicz OC, Bicanic N, Shen FQ. Earthquake input definition and the [65 transmitting boundary conditions. In: St Doltsinis I, editor. Adv. Comput. Nonlinear Mech.. Sydney, Australia: Balkema Publishers; 1988. p. 109-38.

Chung J, Hulbert GM. A time integration algorithm for structural dynamics with improved numerical dissipation: the generalised-a method. J Appl Mech 1993;60:371-5. http://dx.doi.org/10.1115/1.2900803.

Kontoe S. Development of time integration schemes and advanced boundary conditions for dynamic geotechnical analysis PhD thesis. London: Imperial College London; 2006.

Kontoe S, Zdravkovic L, Potts DM. An assessment of time integration schemes for dynamic geotechnical problems. Comput Geotech 2008;35:253-64. http:// dx.doi.org/10.1016/j.compgeo.2007.05.001.

Taborda DMG. Development of constitutive models for application in soil dynamics PhD thesis. London: Imperial College London; 2011. Taborda DMG, Zdravkovic L, Kontoe S, Potts DM. Computational study on the modification of a bounding surface plasticity model for sands. Comput Geotech 2014;59:145-60. http://dx.doi.org/10.1016/j.compgeo.2014.03.005. Papadimitriou AG, Bouckovalas GD. Plasticity model for sand under small and large cyclic strains: a multiaxial formulation. Soil Dyn Earthq Eng 2002;22:191-204. http://dx.doi.org/10.1016/S0267-7261(02)00009-X. Manzari MT, Dafalias YF. A critical state two-surface plasticity model for sands. Géotechnique 1997;47:255-72. http://dx.doi.org/10.1680/geot.1997.47.2.255. Klokidi M. Numerical simulation of liquefaction behaviour with focus on the re-liquefaction of sand deposits MSc thesis. Imperial College London; 2015. Darendeli MB. Development of a new family of normalised modulus reduction and material damping curves PhD thesis. Austin: University of Texas; 2001. GNS. GeoNet 2014 <http://www.geonet.org.nz> [accessed June 15, 2014]. Elgamal A-W, Zeghal M, Tang HT, Stepp JC. Lotung downhole array. I: Evaluation of site dynamic properties. J Geotech Eng 1995;121:350-62. Been K, Jefferies MG. A state parameter for sands. Géotechnique 1985;35:99-112. http://dx.doi.org/10.1680/geot.1985.35.2.99. Wotherspoon LM, Orense R, Bradley BA, Cox BR, Wood C, Green RA. Geotechnical characterization of Christchurch strong motion stations V3. Earthquake Commission Report, Project No. 12/629, University of Canterbury. Civil and Natural Resources Engineering; 2015

Robertson PK. Interpretation of in-situ tests-some insights. In: James K. Mitchell Lect. Proc. 4th Int. Conf. Geotech. Geophys. Site Charact., Porto de Galinhas, Pernambuco, Brazil; 2012, p. 1-22.

NZGD. New Zealand Geotechnical Database - MBIE & EQC 2016 <www.nzgd. org.nz> [accessed September 30, 2016].

Ishihara K. Soil behaviour in earthquake geotechnics. Oxford Engineering Series. Oxford: Oxford University Press; 1996.

Tsai Y-B, Huang M-W. Strong ground motion characteristics of the Chi-Chi, Taiwan, earthquake of September 21, 1999. Earthq Eng Eng Seismol 2000;2:1-21.

Kramer SL. Geotechnical earthquake engineering. New Jersey: Prentice Hall; 1996.

Bozorgnia Y, Campbell KW. The vertical-to-horizontal response spectral ratio and tentative procedures for developing simplified V/H and vertical design spectra. J Earthq Eng 2004;8:175-207.

Papazoglou AJ, Elnashai AS. Analytical and field evidence of the damaging effect of vertical earthquake ground motion. Earthq Eng Struct Dyn 1996;25:1109-38.

Kunnath SK, Erduran E, Chai YH, Yashinsky M. Effect of near-fault vertical ground motions on seismic response of highway overcrossings. J Bridge Eng 2008;13:282-90. http://dx.doi.ore/10.1061/(ASCE)1084-0702(2008)13:3, (282).

Riches LK. Observed earthquake damage to Christchurch city council owned retaining walls and the repair solutions developed. In: Proc. 6th Int. Conf. Earthq. Geotech. Eng., 1-4 November, Christchurch, New Zealand; 2015. Han B. Hydro-mechanical coupling in numerical analysis of geotechnical structures under multi-directional seismic loading PhD thesis. London: Imperial College London; 2014.

Thomas J. Static, cyclic and post liquefaction undrained behaviour of Fraser River Sand MSc thesis. University of British Columbia; 1992. Ghafghazi M. Towards comprehensive interpretation of the state parameter from cone penetration testing in cohesionless soils PhD thesis. University of British Columbia; 2011.

Sriskandakumar S. Cyclic loading response of Fraser River Sand for validation of numerical models simulating centrifuge tests MSc thesis. The University of British Columbia; 2004.

Masing G. Eigenspannungen und verfestigung beim messing. In: Proc. 2nd Int. Conf. Appl. Mech. Zurich, Switz..