High-frequency ground motion ampliﬁcation during the 2011 Tohoku earthquake explained by soil dilatancy

SUMMARY 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 2 g . 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 M w 9.1 earthquake was characterized by increased ampliﬁcation at frequencies above 10 Hz and in peak ground acceleration. This behaviour contrasts with the more common non-linear response encountered at non-liqueﬁable sites, which results in deampliﬁcation at higher frequencies. We simulate propagation of SH waves through the dense sand deposit using a non-linear ﬁnite 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 ﬁnite 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.


I N T RO D U C T I O N
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 in hysteretic 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 ).
However, other observations exist which show that soils do not reduce high-frequency ground motions in all cases. During the 2008 M w 6.9 Iwate-Miyagi earthquake in Japan, the peak acceleration above the fault reached 3.8 × 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 M w 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. Similar dilation pulses which amplitudes between 0.5 × g and 0.7 × g were recorded at the stations NNBS, REHS and CBGS during the 2011 M w 6.3 Christchurch earthquake (e.g. Bradley & Cubrinovski 2011). During the 2011 M w 9.1 Tohoku earthquake, which caused widespread liquefaction, dilation pulses were also widely observed, for example, at the K-NET stations CHB024 and MYG013 (Bonilla et al., 2011b), with a dilation pulse exceeding 1 × 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 × g (Fig. 1).

O P R E C O R D S O F T H E 0 1 1 T O H O K U E A RT H Q UA K E
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.5 g) 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.13 m s −2 (2 × 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 M w 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 (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.

WAV E P RO PA G AT I O N S I M U L AT I O N S
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 Average response spectral ratios calculated from 11 events with M w 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 M w 9.1 event. Table 1. Onahama port soil model.

Top
Description vector. We chose a spatial discretization of 0.10 m and a time step of 1.05 × 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 accelerometer is located inside a silt layer that begins 7 m below the surface (Table 1). Shearvelocities v s and densities ρ were adopted from Yamazaki & Gotoh (2011). The coefficient of Earth at rest K 0 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 N-values based on empirical relations (Terzaghi et al. 1996).

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 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 σ m versus shear stress τ xy ) is controlled by the liquefaction front parameter S 0 , which is a function of the cumulative shear work w: The parameter w 1 defines the contribution of normalized shear work over the entire zone of S 0 , while p 1 controls the initial phase of dilatancy and p 2 the final phase. The parameter S 1 , typically set to 0.01, is required for numerical stability and prevents S 0 from becoming zero, that is, liquefaction is not allowed. Additionally, the threshold ratio c 1 is required, which describes the threshold of shear work for which no pore water pressure buildup occurs. The dilatancy parameters w 1 , p 1 , p 2 and c 1 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: where, SA (k) o and SA (k) s represent the observed and simulated spectral accelerations at frequency k, while U (n) o and U (n) s indicate the observed and simulated upper (positive) envelopes at time step n; conversely L (n) o and L (n) s 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 U (n) 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 w 1 and the threshold ratio c 1 in both sand layers (layer 2 and layer 3). Since the misfit is not very sensitive to the choice of p 1 and p 2 , we constrained p 1 to 0.6 and p 2 to 1.2 in the liquefiable sands (Table 1). To further reduce the number of unknown parameters we defined w 1 = 1000 for the lower sand layer. This value was chosen based on the SPT N-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 w 1 may vary over several orders of magnitude, we inverted for log 10 c and log 10 w 1 . 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.

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 w 1 (Fig. 3b) are well constrained. Low misfits (<0.01) are obtained for log 10 c > 4 and 1.5 < log 10 w 1 < 2.5. The misfit is less sensitive to the choice of c 1 in the upper sand layer (layer 2, Fig. 3c), but lower misfits are generally obtained for smaller values of c 1 < 2.5. In the lower sand layer (layer 3, Fig. 3d), low misfits are obtained for values of c 1 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 × 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 U max as a function of depth for all the models sampled during the inversion. Inside the upper sand layer U max approaches the initial effective mean stress σ 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 U max remains well below σ 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 & Elgamal, 1994;Holzer & Youd, 2007), which are recorded as dilation pulses on the surface.

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 M w 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 w 1 ≈ 115 obtained for the upper sand layer by inversion is much higher than the values determined for other sites, where w 1 typically ranges between 1.5 and 7 (e.g. Iai et al. 1990bIai et al. , 1995Roten 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.

C O N C L U S I O N S
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 × g have been recorded on soft soils which responded distinctively non-linear to the shaking. Such observations suggest that careful consideration should be given to pore water pressure effects when evaluating the role of soil nonlinearity 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 highfrequency 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.

A C K N O W L E D G E M E N T S
Strong motion records were provided by the Port and Airport Research Institute (Nozu, 2004; http://www.eq.pari.go.jp/ 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 http://www.openstreetmap.org ( C OpenStreetMap contributors, CC BY-SA). This research is partly funded through a contract with the Swiss Nuclear Safety Inspectorate (ENSI).