Scholarly article on topic 'High-frequency ground motion amplification during the 2011 Tohoku earthquake explained by soil dilatancy'

High-frequency ground motion amplification during the 2011 Tohoku earthquake explained by soil dilatancy Academic research paper on "Earth and related environmental sciences"

Share paper

Academic research paper on topic "High-frequency ground motion amplification during the 2011 Tohoku earthquake explained by soil dilatancy"

Geophysical Journal International

Geophys. J. Int. (2013) 193, 898-904

doi: 10.1093/gji/ggt001

High-frequency ground motion amplification during the 2011 Tohoku earthquake explained by soil dilatancy

D. Roten,1 D. Fah1 and L. F. Bonilla2

1 Swiss Seismological Service, Sonneggstrasse 5, 8092 Ziirich, Switzerland. E-mail:

2 Université Paris-Est, IFSTTAR, 58 Bd Lefebvre 75732 Paris Cedex 15, France

Accepted 2013 January 2. Received 2012 September 25; in original form 2012 June 21

'au IS)


Ground motions of the 2011 Tohoku earthquake recorded at Onahama port (Iwaki, Fukushima prefecture) rank among the highest accelerations ever observed, with the peak amplitude of the 3-D acceleration vector approaching 2g. The response of the site was distinctively non-linear, as indicated by the presence of horizontal acceleration spikes which have been linked to cyclic mobility during similar observations. Compared to records of weak ground motions, the response of the site during the Mw 9.1 earthquake was characterized by increased amplification at frequencies above 10 Hz and in peak ground acceleration. This behaviour contrasts with the more common non-linear response encountered at non-liquefiable sites, which results in deamplification at higher frequencies. We simulate propagation of SH waves through the dense sand deposit using a non-linear finite difference code that is capable of modelling the development of excess pore water pressure. Dynamic soil parameters are calibrated using a direct search method that minimizes the difference between observed and simulated acceleration envelopes and response spectra. The finite difference simulations yield surface acceleration time-series that are consistent with the observations in shape and amplitude, pointing towards soil dilatancy as a likely explanation for the high-frequency pulses recorded at Onahama port. The simulations also suggest that the occurrence of high-frequency spikes coincided with a rapid increase in pore water pressure in the upper part of the sand deposit between 145 and 170 s. This sudden increase is possibly linked to a burst of high-frequency energy from a large slip patch below the Iwaki region.

Key words: Earthquake ground motions; Site effects; Computational seismology.


Understanding how near-surface geological conditions affect ground motions at a given site is crucial for both interpreting acceleration time-series and for seismic hazard assessment. The linear response of soils to weak ground motions, which consists in amplification of seismic waves and increased shaking duration, is well established and can be attributed to the impedance contrast between the soft soil and the underlying bedrock (e.g. Kramer 1996). However, laboratory testing shows that soils exhibit a decrease in shear stiffness and an increase inhysteretic damping when sheared to large strains (e.g. Seed & Idriss 1969). Field evidence of this non-linear behaviour is represented by observations of reduced amplification during strong ground motions compared to weak ground motions, especially at high frequencies (f > 1 Hz), and by a shift in the resonance frequency of the soil deposit (e.g. Beresnev & Wen 1996). Examples of such observations include the 2000 M 7.3 Tottori earthquake recorded at the site TTRH02 (Bonilla et al. 2011a), or the 2011 M 9.1 Tohoku earthquake recorded at the stations IBRH11 and IBRH16 (Regnier et al. 2011).

However, other observations exist which show that soils do not reduce high-frequency ground motions in all cases. During the 2008 Mw 6.9 Iwate-Miyagi earthquake in Japan, the peak acceleration above the fault reached 3.8 x g in the upward direction, with more than twice the acceleration recorded in the upward direction than in the downward direction. Aoi et al. (2008) explained these observations with a trampoline-like effect in loose soils, with the highest acceleration caused by rebound of particles following a quasi freefall state. Similar asymmetric waveforms were observed during the 2011 Mw 6.3 Christchurch earthquake (e.g. Fry et al. 2011).

In the horizontal directions, high-frequency acceleration spikes can be induced by cyclic mobility (i.e. buildup of pore-water pressure) in soils susceptible to liquefaction (Iai et al. 1995; Archuleta et al. 2000; Bonilla et al. 2005). The observations made at the Wildlife liquefaction array during the 1986 Superstition hills earthquake are a well-known example for this behaviour (Holzer et al. 1989). Zeghal & Elgamal (1994) linked the high-frequency spikes in the Wildlife acceleration data to temporary drops in the recorded pore water pressure, which are caused by the dilative behaviour of sands and result in a temporary recovery of shear strength.

© The Authors 2013. Published by Oxford University Press on behalf of The Royal Astronomical Society.


-2.0 m -s"2

2.4 m-s"

Figure 1. Acceleration time-series recorded at Onahama port (Fig. 2) in the borehole at 11 m depth (GL-11 m) and on the surface (GL-0 m). Numbers above the traces indicate peak ground acceleration (PGA).

Similar dilation pulses which amplitudes between 0.5 x g and 0.7 x g were recorded at the stations NNBS, REHS and CBGS during the 2011 Mw 6.3 Christchurch earthquake (e.g. Bradley & Cubrinovski 2011). During the 2011 Mw 9.1 Tohoku earthquake, which caused widespread liquefaction, dilation pulses were also widely observed, for example, at the K-NET stations CHB024 andMYG013 (Bonilla et al., 2011b), with a dilation pulse exceeding 1 x g at the latter site. In this paper we analyse the Onahama port (OP) records of the Tohoku earthquake (Wakai & Nozu 2011), which are characterized by high-frequency acceleration spikes exceeding 1.5 x g (Fig. 1).

(100 Hz). During the 2011 Tohoku earthquake, however, the soil deposit amplified ground motions by a factor of 9.5 at 28 Hz and by a factor of 7 at 100 Hz (Fig. 2b). This increase in high-frequency amplification is the opposite of the non-linear response that is typically associated with non-liquefiable geomaterials, such as clays or silts, which results in a reduction of high-frequency energy due to hysteretic damping. In the next section we perform numerical simulations which show that the OP records may be reproduced using an advanced constitutive soil model that is capable of predicting pore pressure fluctuations.


The acceleration time-series shown in Fig. 1 were recorded at OP near Iwaki (Fukushima prefecture), on a vertical array operated by the Port and Airport Research Institute (Fig. 2a). The borehole accelerometer 11 m below the surface recorded a peak ground acceleration (PGA) of 1.7 m s-2 in the E-W and 2.0 m s-2 in the N-S direction (Fig. 1). Acceleration time-series on the surface are characterized by high-frequency pulses, which result in PGAs of 15.3 m s-2 (1.5g) in the E-W direction and 11.4 m s-2 (1.1g) in the N-S direction. The maximum amplitude of the 3-D acceleration vector is 19.13ms-2 (2 x g), which ranks among the highest values ever recorded (Anderson, 2010).

We computed average response spectral ratios between the surface and borehole receiver for 11 Mw 4.5-7.6 earthquakes recorded at OP (Fig. 2a). These events indicate an average amplification of ~5.7 at 5.4 Hz (Fig. 2b) and an amplification of ~4 in PGA


Using the fully non-linear 1-D finite difference (FD) code NOAH (Bonilla et al. 2005) we simulate the propagation of SH waves through 11 m of horizontally layered soil. Simulations are performed using borehole boundary conditions, which means that the ground motion at the base of the soil column is prescribed by the recorded borehole signal at any time (Joyner & Chen 1975). This type of boundary condition allows to use the downhole record, which contains both the incident and reflected wavefield, directly as input ground motion. Elastic (transmitting) boundary conditions, on the other hand, would require an input signal that contains the incident wavefield only. Therefore, borehole boundary conditions are often used for non-linear site response studies (e.g. Iai et al. 1995; Bonilla et al. 2005), because it is not possible to separate the incoming from the reflected wavefield after strong non-linearity has taken place. We perform simulations for horizontal ground motions in the direction N126°E to capture the peak in the acceleration

Figure 2. Map of eastern Honshu with epicentres of the events used for calculation of response spectral ratios. The red rectangle (south of Iwaki) shows the extent of the inset map. The location of the vertical array operated by the Port and Airport Research Institute is indicated by the yellow triangle. (b) Average response spectral ratios calculated from 11 events with Mw ranging from 4.5 to 7.6 (thick black line) ±standard deviation (grey shaded area). The red line shows the response spectral ratio for the Mw 9.1 event.

Table 1. Onahama port soil model.

Top Description Low-strain Non-linear

(m) parameters parameters

0 Fill soil vs = 100 ms-1 V = 30°

unsaturated p = 1800 kg m-3 £ max = 20 percent

(layer 1) K0 = 1.0 c

i.25 Upper sand v s = 124 m s- 1 V = 30°

(layer 2) p = 2000 kg m-3 Vp = 20°

K0 = 1.0 W1, C1

n = 40 per cent P1 = 0.6, p2 = 1.2

3.50 Lower sand vs = 215 ms-1 V = 40°

(layer 3) p = 2000 kg m-3 Vp = 28°

K0 = 1.0 w1 = 1000, c1

n = 40 per cent p1 = 0.6, p2 = 1.2

7 Silt vs = 950 ms-1 V = 20°

(layer 4) p = 2200 kg m-3 £ max = 20 percent

K0 = 1.0 c=103 N m-2

vector. We chose a spatial discretization of 0.10 m and a time step of 1.05 x 10-5 s, which allows us to model frequencies up to 20 Hz using 50 gridpoints per wavelength.

The deposits at the OP site consist of 1.25 m of fill soil located above the ground water table (Yamazaki & Gotoh 2011). Below the fill soil two saturated sand layers of 2.25 and 3.50 m thickness, respectively, are encountered. The borehole accelerometeris located inside a silt layer that begins 7 m below the surface (Table 1). Shear-velocities vs and densities p were adopted from Yamazaki & Gotoh (2011). The coefficient of Earth at rest K0 was assumed to be equal 1 for all layers (Table 1). We determined the friction angles cohesion c and the maximum damping ratio §max from the reported SPT ^-values based on empirical relations (Terzaghi et al. 1996).

mobility and soil dilatancy (Iai et al. 1990a). The key concept behind the Iai model is the liquefaction front, which describes the decrease of effective mean stress due to the increase in pore water pressure. The position of the liquefaction front in normalized stress space (effective confining stress am versus shear stress Txy) is controlled by the liquefaction front parameter S0, which is a function of the cumulative shear work w:

if (w < wi),

S0 = i - 0.6 ( — wi

So = (0.4 - Si) — + Si if(w>wi).

The parameter w1 defines the contribution of normalized shear work over the entire zone of S0, while p1 controls the initial phase of dilatancy andp2 the final phase. The parameter S1, typically set to 0.01, is required for numerical stability and prevents S0 from becoming zero, that is, liquefaction is not allowed. Additionally, the threshold ratio c1 is required, which describes the threshold of shear work for which no pore water pressure buildup occurs. The dilatancy parameters w1, p1, p2 and c1 are typically calibrated from stress-controlled laboratory experiments (e.g. Iai et al. 1990b; Roten et al. 2009). For this study, however, we defined the dilatancy parameters directly from the strong motion records.

We performed a direct search using the Neighbourhood algorithm (Sambridge, 1999) to find a set of dilatancy parameters that reproduces the observed horizontal accelerations at the surface (Fig. 1). The misfit is defined as the least-square difference between the observed and simulated response spectra and the observed and simulated envelopes:

i A (SAo - SA )2

misfit = —

|SAo |

(L0n) - LSn))2

3.1 Calibration of dilatancy parameters

NOAH applies effective stress analysis for liquefiable layers based on the multispring mechanism model introduced by Towhata & Ishihara (1985), and an extension of this model which treats cyclic

U0n) - Un))

where, SAOk) and SA® represent the observed and simulated spectral accelerations at frequency k, while U0n) and USn) indicate the

Table 2. Inverted soil parameters and minimum misfit value.

Parameter Layer Symbol Unit Range Value

Cohesion 1 logioc Nm-2 0-7 6.04

Overall dilatancy 2 logiowi - -1-5 2.06

Threshold limit 2 ci - 0.01-5.0 0.32

Threshold limit 3 ci - 0.01-20.0 16.23

observed and simulated upper (positive) envelopes at time step n; conversely L0n) and L(n) indicate the lower (negative) envelopes. We introduced envelopes to make the misfit more tolerant against small time differences between simulated and observed dilation pulses. The upper and lower envelopes Un) and L(n) at time step n were obtained by multiplying the signal with a Hanning window of 1.5 s width centred at n and determining the maximum and minimum of the resulting time-series.

We inverted for the overall dilatancy w1 and the threshold ratio c1 in both sand layers (layer 2 and layer 3). Since the misfit is not very sensitive to the choice ofp1 andp2, we constrainedp\ to 0.6 and p2 to 1.2 in the liquefiable sands (Table 1). To further reduce the number of unknown parameters we defined w1 = 1000 for the lower sand layer. This value was chosen based on the SPT ^-values above 50 and the high fines content (FC > 15 per cent) reported for this unit (Yamazaki & Gotoh, 2011), which precludes a significant buildup of pore water pressure.

However, we found that the misfit is sensitive to the strength of the uppermost layer, as the amplitudes of dilation pulses originating in the liquefiable sand are reduced by hysteretic damping inside the fill soil before they reach the surface. Therefore, we introduced a cohesion c for the fill soil (layer 1). Because c and w1 may vary over several orders of magnitude, we inverted for log10c and logi0wi. Inverted parameters and the sampled range are listed in Table 2. We performed 30 iterations using a sample size of 24 models, resampling 12 models during each iteration.

3.2 Inversion results

Figs 3(a)-(d) show the model misfit as a function of parameter value for the inversion of the N126°E acceleration time-series. The cohesion c of the fill soil (Fig. 3a) and the overall dilatancy in the upper sand layer w1 (Fig. 3b) are well constrained. Low misfits (<0.01) are obtained for log10c > 4 and 1.5 < log10w1 < 2.5. The misfit is less sensitive to the choice of c1 in the upper sand layer (layer 2, Fig. 3c), but lower misfits are generally obtained for smaller values of c1 < 2.5. In the lower sand layer (layer 3, Fig. 3d), low misfits are obtained for values of ci around 2 and 16, indicating that this parameter has little influence on the overall misfit. The values of the minimum misfit model are given in Table 2. Envelopes of the simulated acceleration time-series for the minimum misfit model (Fig. 3e) show similar features as the observed envelopes. The observed response spectra (Fig. 3f) is also well reproduced by the minimum misfit model. The strongest pulse in the simulation reaches ~15.9 m s-2 (1.6 x g) at ~157.6 s (relative to hypocentral time), at the same time as in the observed time-series (Fig. 4b).

Fig. 5(a) shows the maximum excess pore water pressure Umax as a function of depth for all the models sampled during the inversion. Inside the upper sand layer Umax approaches the initial effective mean stress a'm0 (dashed line in Fig. 5a) for all models with misfits below ~0.01. Pore water pressure buildup is also occurring in the lower sand layer, but Umax remains well below a'm0 during the

shaking. Models with misfits below ~0.01 indicate that buildup of excess pore water pressure in the upper sand layer occurred mostly between 145 and 170 s (Fig. 5b), which coincides with the time when high-amplitude dilation pulses were recorded on the surface. Models which result in a faster buildup of pore water pressure are characterized by a higher misfit, as well as models which produce no significant buildup of pore water pressure. The stress path of the minimum misfit model (Fig. 5c) repeatedly crosses the phase transformation line between 145 and 180 s, which indicates that the soil enters dilative behaviour. In this state the stress-strain loop (Fig. 5d) follows an S-shape and develops large strains up to 4 per cent, followed by large and spiky shear stresses (e.g. Zeghal & El-gamal, 1994; Holzer & Youd, 2007), which are recorded as dilation pulses on the surface.

3.3 Discussion

The generally good agreement between the minimum misfit model and the observations suggests that soil dilatancy is a likely explanation for the high-frequency pulses recorded at OP. Even though no pore water pressure transducers are available for the OP site, the inversion for dilatancy parameters (Fig. 5b) reveals a possible pattern of pore water pressure development that would explain the data. Strong motion stations in the Ibaraki region recorded a burst of high-frequency energy between 130 and 170 s, which has been attributed to a large slip patch beneath the northern end of Ibaraki either over the plate interface or in the crust (e.g. Ide et al. 2011; Furumura et al. 2011). Such a slip patch has also been identified in various back-projection studies (e.g. Meng et al. 2011; Roten et al. 2012). This strongly suggests that the rapid increase in pore water pressure between 145 and 170 s (Fig. 5b), which resulted in dilation pulsed during the same time range, was triggered by high-frequency energy originating from this asperity below the Ibaraki region.

While significant buildup of pore water pressure is required for explaining the high-frequency pulses with dilatancy, the recorded surface accelerations do not show a lack of high-frequency towards the end of the record. This may indicate that the buildup of pore water pressure has not been sufficient to trigger complete liquefaction of the upper sand layer. A number of sand boils were identified at OP when the station was visited in May 2011 (Atsushi Wakai, personal communication, 2012). However, it is not clear whether these sandboils developed during the main shock or during a powerful aftershock. The tsunami that flooded the area shortly after the Mw 9.1 rupture may have removed sandboils associated with the main event. However, liquefaction-induced damage, such as uneven settlement or tilting, has been reported for numerous structures in Iwaki and northern Ibaraki (Aydan et al. 2011).

The overall dilatancy of w1 & 115 obtained for the upper sand layer by inversion is much higher than the values determined for other sites, where w1 typically ranges between 1.5 and 7 (e.g. Iai et al. 1990b, 1995; Roten et al. 2009). This indicates that the liquefaction resistance of this layer is high, that is, a large amount of shear work is required to build up pore water pressure. The unusually long duration (> 3 min) of the Tohoku earthquake might have provided just enough loading cycles to initiate cyclic mobility. Indeed we did not discover dilation pulses in any other acceleration time-series recorded by this station, and the Tohoku earthquake represents the only event where PGAs exceeded 1.5 m s-2 at OP. Unfortunately the record is incomplete, since the OP station went out of operation a few hours after it was submerged by the tsunami, and many aftershocks that might have triggered cyclic mobility were not recorded.

Figure 3. Inversion of N126° E acceleration time-series with model misfit as a function of parameter value for (a) log10c in layer 1, (b) log10W1 in layer 2, (c) c1 in layer 2, and (d) c1 in layer 3. (e) Simulated envelopes of best-fit model versus observed envelopes; (f) simulated spectral acceleration (SA) of best-fit model versus observed SA.

'in É

observed — minimum misfit model

15 10 5 0 -5 150

152 154

156 158 160

80 100 120 140 160 180

Time (s)

Figure 4. (a) Simulated (dark blue) and observed (grey) surface acceleration time-series. The rectangle shows the enlarged sections given in (b).


This rather extreme example of cyclic mobility illustrates the complexity of the processes which may dramatically alter the response of liquefiable soils during strong shaking. It shows that soil non-

linearity is not always manifest as a reduction in high-frequency amplification, but may also lead to increased amplification at high frequencies and in PGAs in the horizontal directions. As a result accelerations in excess of 1 x g have been recorded on soft soils which responded distinctively non-linear to the shaking. Such

Figure 5. Excess pore water pressure development of models sampled during inversion. (a) Maximum excess water pressure Umax as a function of depth. The colour code reflects misfit, with lower misfit models plotted on top. (b) Cumulative maximum excess pore water pressure as a function of time at 2 m depth in the upper sand layer. Pore water pressure drops caused by dilation are not shown for clarity. Dashed lines show initial effective mean stress Omq. (c) Stress path (including failure line, dashed, and phase transformation line, dotted) and (d) stress-strain loop of best-fit model at 2 m depth.

observations suggest that careful consideration should be given to pore water pressure effects when evaluating the role of soil non-linearity as a physical limit to ground motions (e.g. Frankel, 1999; Hanks et al., 2005) at high frequencies. Even though high-amplitude dilation pulses are rarely observed, they might be relevant for the hazard curve at low likelihoods of exceedance, which is required for critical facilities such as nuclear power stations. For such installations, frequencies that are too high to cause structural damage (>15 Hz) require attention, because the response of essential systems and equipment to high-frequency seismic input needs to be evaluated.

It is obvious that site response analysis neglecting pore water pressure buildup fails to predict the response of sites like OP during strong shaking (e.g. Kramer et al., 2011). The good match between the synthetic and observed surface acceleration presented in this study demonstrates that advanced constitutive models, which treat pore pressure fluctuations, have the capability to predict high-frequency dilation pulses. Calibration of soil models remains a challenge, however, as the dilatancy parameters required in the Iai et al. (1990a) and similar cyclic mobility models have only been determined for very few sites. Direct inversion of strong motion records

for dilatancy parameters promises to identify models that are able to explain observations, and will help to assess the reliability of existing laboratory and field tests to characterize the potential of soils to generate pore water pressure. This highlights the necessity of acquiring more strong motion records using borehole arrays, which should be equipped with pore pressure transducers at different depths.


Strong motion records were provided by the Port and Airport Research Institute (Nozu, 2004; kyosin/, in Japanese). We thank Atsushi Wakai for providing pictures of sandboils at the OP site and for helpful information. Malcolm Sambridge contributed the Neighbourhood algorithm. Fig. 2(a) was generated with GMT (Wessel & Smith, 1998) using map data from (©OpenStreetMap contributors, CC BY-SA). This research is partly funded through a contract with the Swiss Nuclear Safety Inspectorate (ENSI).


Anderson, J., 2010. Source and site characteristics of earthquakes that have caused exceptional ground accelerations and velocities, Bull. seism. Soc. Am., 100(1), 1-36.

Aoi, S., Kunugi, T. & Fujiwara, H., 2008. Trampoline effect in extreme ground motion, Science, 322(5902), 727.

Archuleta, R., Bonilla, L. & Lavallee, D., 2000. Nonlinearity in observed and computed accelerograms, in Paper No. 1934, Proceedings of the 12th World Conference on Earthquake Engineering.

Aydan, O., Iwatate, T. & Tano, H., 2011. Liquefaction-induced damage caused by M 9.0 East Japan mega earthquake on March 11, 2011, JSCE report, Japan Society of Civil Engineers.

Beresnev, I.A. & Wen, K.-L., 1996. Nonlinear soil response—a reality?, Bull. seism. Soc. Am., 86(6), 1964-1978.

Bonilla, L., Archuleta, R. & Lavallee, D., 2005. Hysteretic and dilatant behavior of cohesionless soils and their effects on nonlinear site response: field data observations and modeling, Bull. seism. Soc. Am., 95(6), 23732395.

Bonilla, L., Gelis, C. & Regnier, J., 2011a. The challenge of nonlinear site response: field data observations and numerical simulations, in 4th IASPEI/IAEE International Symposium: Effects of Surface Geology on Seismic Motion, August 23-26, 2011, University of Santa Barbara.

Bonilla, L., Tsuda, K., Pulido, N., Regnier, J. & Laurendeau, A., 2011b. Nonlinear site response evidence of K-NET and KiK-net records from the 2011 off the Pacific coast of Tohoku Earthquake, Earth Planets Space, 63(7), 785-789.

Bradley, B. & Cubrinovski, M., 2011. Near-source strong ground motions observed in the 22 February 2011 Christchurch earthquake, Seism. Res. Lett., 82(6), 853-865.

Frankel, A., 1999. How does the ground shake?, Science, 283(5410), 20322033.

Fry, B., Benites, R. & Kaiser, A., 2011. The Character of accelerations in the Mw 6.2 christchurch earthquake, Seism. Res. Lett., 82(6), 846-852.

Furumura, T., Takemura, S., Noguchi, S., Takemoto, T., Maeda, T., Iwai, K. & Padhy, S., 2011. Strong ground motions from the 2011 off-the Pacific-Coast-of-Tohoku, Japan (Mw = 9.0) earthquake obtained from a dense nationwide seismic network, Landslides, 333-338, 1-6.

Hanks, T., Abrahamson, N., Board, M., Boore, D., Brune, J. & Cornell, C., 2005. Observed ground motions, extreme ground motions, and physical limits to ground motions, in Directions in Strong Motion Instrumentation, Polat Giilkan, John G. Anderson eds, pp.

55-59. Springer-Verlag, Berlin.

Holzer, T. & Youd, T., 2007. Liquefaction, ground oscillation, and soil deformation at the Wildlife Array, California, Bull. seism. Soc. Am., 97(3), 961-976.

Holzer, T., Youd, T.L. & Hanks, T., 1989. Dynamics of Liquefaction during the 1987 Superstition Hills, California, Earthquake, Science, 244(4900),


Iai, S., Matsunaga, Y. & Kameoka, T., 1990a. Strain space plasticity model for cyclic mobility, Report Port Harb. Res. Inst., 29, 27-56.

Iai, S., Matsunaga, Y. & Kameoka, T., 1990b. Parameter identification for a cyclic mobility model, Report Port Harb. Res. Inst., 29, 57-83.

Iai, S., Morita, T., Kameoka, T., Matsunaga, Y. & Abiko, K., 1995. Response of a dense sand deposit during 1993 Kushiro-Oki Earthquake, Soils Found., 35, 115-131.

Ide, S., Baltay, A. & Beroza, G., 2011. Shallow dynamic overshoot and energetic deep rupture in the 2011 Mw 9.0 Tohoku-Oki earthquake, Science, 332(6036), 1426-1429.

Joyner, W. & Chen, A., 1975. Calculation of nonlinear ground response in earthquakes, Bull. seism. Soc. Am., 65(5), 1315-1336.

Kramer, S.L., 1996. Geotechnical earthquake engineering, Prentice Hall, New Jersey.

Kramer, S.L., Hartvigsen, A.J., Sideras, S.S. & Ozener, P.T., 2011. Site response modeling in liquefiable soil deposits, in 4th IASPEI/IAEE International Symposium: Effects of Surface Geology on Seismic Motion, August 23-26, 2011, University of Santa Barbara.

Meng, L., Inbal, A. & Ampuero, J., 2011. A window into the complexity of the dynamic rupture of the 2011 Mw 9 Tohoku-Oki earthquake, Geophys. Res. Lett., 38, L00G07, doi:10.1029/2011GL048118.

Nozu, A., 2004. Current status of strong-motion earthquake observation in Japanese ports, Journal of Japan Association of Earthquake Engineering, 4(3), 79-83.

Regnier, J., Bonilla, L., Betrand, E. & Semblat, J., 2011. Empirical evidence of nonlinear site response at several KiK-net stations, in 4th IASPEI/IAEE International Symposium: Effects of Surface Geology on Seismic Motion, August 23-26, 2011, University of Santa Barbara.

Roten, D., Fah, D., Bonilla, L., Alvarez-Rubio, S., Weber, T. & Laue, J., 2009. Estimation of non-linear site response in a deep Alpine valley, Geophys. J. Int., 178(3), 1597-1613.

Roten, D., Miyake, H. & Koketsu, K., 2012. A Rayleigh wave back-projection method applied to the 2011 Tohoku earthquake, Geophys. Res. Lett., 39(2), L02302, doi:10.1029/2011GL050183.

Sambridge, M., 1999. Geophysical inversion with a neighbourhood algorithm: I. searching a parameter space, Geophys. J. Int., 138, 479494.

Seed, H. & Idriss, I., 1969. Influence of soil conditions on ground motions during earthquakes, J. SoilMech. Found. Div., ASCE, 95(SM1), 99137.

Terzaghi, K., Peck, R. & Mesri, G., 1996. Soil Mechanics in Engineering Practice, Wiley-Interscience.

Towhata, I. & Ishihara, K., 1985. Modeling soil behavior under principal axes rotation, in Fifth International Conference on Numerical Methods in Geomechanics, Nagoya, Japan, pp. 523-530.

Wakai, A. & Nozu, A., 2011. The strong motion observation network in Japanese ports, in NEES workshop: Using the NEES Equipment Site Facilities in ESG Research and International Collaborations, August 23, 2011, University of Santa Barbara.

Wessel, P. & Smith, W., 1998. New, improved version of generic mapping tools released, EOS, Trans. Am. geophys. Un., 79(47), 579.

Yamazaki, H. & Gotoh, Y., 2011. Study on liquefaction of soil subjected to seismic motion in the 2011 off the Pacific coast of Tohoku earthquake, Technical note no. 1242, The Port and Airport Research Institute, (in Japanese with English abstract).

Zeghal, M. & Elgamal, A., 1994. Analysis of site liquefaction using earthquake records, J. Geotech. Eng.;(United States), 120(6), 996-1017.