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

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/). ⇑ 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). Vasiliki Tsaparli a,⇑, Stavroula Kontoe , David M.G. Taborda , David M. Potts a


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][2][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][8][9][10][11][12][13] and, as such, more EBMs for liquefaction triggering have been developed [10,11,[13][14][15][16][17][18][19][20][21]. 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 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 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.

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 Table1.
The maximum shear modulus, G max , is given by Eq. (1) [33]: where B is the elastic shear modulus constant and p 0 ref a reference pressure. Assuming a constant Poisson's ratio, the average compressional 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/s 2 , 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 hydromechanically 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 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 compressional 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 Â 0.1 m 2 . 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, q 1 , 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 timehistories and Fourier Spectra for each of the two motions considered, while Fig. 2 depicts the corresponding initial and final mean  (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.

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.

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 Christchurch 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 M w 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 timehistories 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 (f 1 ) 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 f 1 , 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 0 -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 postresonance 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 energy time-histories have been based on vertical effective stressvertical strain loops and have been calculated using the algorithm implemented in ICFEP in three-dimensional stress space by Taborda et al. [27]. The simplification for one-dimensional loading conditions of the general methodology for the calculation of continuous histories of energy dissipation and damping is shown    schematically in Fig. 8 and is applicable to irregular motions where the stress-strain loops are not necessarily closed.
The accumulated energy, E acc;i , 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): 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): 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: 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 f 1 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, r p , has been defined in Tsaparli et al. [28] according to Eq. (7): where Dp 0 is the change in mean effective stress since the start of the analysis and p 0 0 is the mean effective stress level after consolidation, prior to the application of the dynamic loading. r p replaces the commonly used excess pore water pressure ratio, r u , defined as the ratio of the change in pore water pressure, Du, over the vertical effective stress after consolidation, r 0 v0 , 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 r u values greater than unity even when liquefaction has not taken place. Liquefaction occurrence is defined here by an r p value equal to about 0.9, with a value of unity corresponding to complete loss of soil strength (initial liquefaction [54]), similar to r u . Fig. 10 shows r p variations against time (t), dissipated energy (E diss ), number of response cycles (N tot ), representative cycles based on stress ðN stress;repr Þ and representative cycles based on energy (N energy;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, E diss;loop and E el;loop , are calculated and stored as described earlier (Eqs. (4) and (5)) and schematically  shown in Fig. 8, whenever two consecutive reversals, rev i and rev iþ1 , are registered. As such, the hysteretic damping ratio of a cycle is given by Eq. (8): 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, E diss;loop , larger than or equal to 30% of the maximum stress amplitude or E max diss;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 r p increases from the early stages of seismic loading in LV_II_166, reaching values of 1 towards the end of the 20 s duration, which indicates the occurrence of liquefaction. Conversely, the value of r p 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 r p time-histories, due to the total stress changes and the assumption of compressible pore fluid. The maximum value that r p 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 r p is plotted against dissipated energy, E diss , it becomes clear that the results from both analyses initially coincide, however, as anticipated, more energy is dissipated overall in the case of LV_II_166. The representation of r p against the response cycles, N tot , 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 0 , takes place. Finally, when r p is plotted against the number of representative cycles based on stress, N stress;repr , or energy, N energy;repr , the variations from the two analyses once again initially are very similar. Analysis LV_II_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.

Number of cycles of the input motion
It 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 M w = 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 ICFEP, 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. It 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, a max , whereas the Chi-Chi one is equivalent to 18.5 cycles with an amplitude of 65% of a max, 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.  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.

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. It 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 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    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 timehistories described earlier [27]. Two different values of damping ratio have been calculated: The first one, denoted as D 1 , is the effective stress mechanism based on vertical effective stress-vertical strain (r 0 v À e v Þ loops; this reflects the solid skeleton. The second one, denoted as D 2 , is the total stress mechanism based on vertical total stress-vertical strain (r v À e v Þ loops; this accounts for the presence of the non-hysteretic water in the soil's pores. 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 (D 1 ) 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 D 1 after about 10 s shown in Fig. 15. Conversely, the computed damping ratio based on total stress components (D 2 ) retains significantly low values justifying the amplification seen towards the ground surface. Average values of D 2 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 D 2 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 compressional structural damage, as this has been observed in a number of field observations [57][58][59][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 Â 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.

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 postresonance 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 Institute 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. It is shown that the model, through the state parameter concept, is able to capture successfully the mechanical behaviour of FRS under different initial states. It 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  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 D r 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, c c , 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 (p 0 0 = 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.