Abstract
Liquefaction-induced lateral spreading during Tohoku earthquake resulted in significant damage, and disruption of functionality for structures and life. The paper aims at reproducing this on-site evidence presenting the state of the art about the most credited qualitative approach and comparing these methods with numerical computation. In this regard, finite element (FE) simulations are increasingly providing a versatile environment in order to assess economical and effective damage. In the study, several systematic three-dimensional FE computations have been conducted to numerically evaluate the effects in terms of liquefaction-induced lateral deformations. The analysis is performed in correspondence with Urayasu City, where the registered liquefaction consequences on residential buildings were wide if compared with the ordinary seismic shake. This study can be used both for post-earthquake evaluations and for pre-earthquake vulnerability predictions.
1. Introduction
The 11th of March 2011 (05:46:24.51 UTC) Tohoku Earthquake can be considered the largest ( 9.0-9.1) in the recent history of Japan and one of the five largest earthquakes of the modern era (United States Geological Survey’s USGS). The earthquake excited a large tsunami that devastated coastal communities in Japan. The most sites affected by liquefaction-induced deformations are inside the named Kanto Plate, a very recent zone with alluvial formations on which many Japanese cities stand. The main effects of liquefaction-induced ground failures were observed around the northern and north-eastern shorelines of Tokyo Bay (e.g., Shin Kiba, Urayasu, Inage, Kaihin Makuhari, Chiba, Isobe, and Mihama), communities along the river Tone (Choshi, Sawara, Itako, Katori, and Kamisu), and areas along the Naka River including Hitachinaka, Miko, and Oarai [1–6].
In particular, the main damages were observed in the Kanto Plain region at many urban cities, which includes the Tokyo Bay and Tone River areas. The liquefied soils were fill materials or young alluvium [5]. Lateral spreading consisted of the development of large horizontal ground displacements due to earthquake-induced liquefaction. This phenomenon resulted in significant damage and considerable replacement costs for existing buildings and civil engineering structures (quay walls, bridge piers, etc.) since it imposed notable lateral loads and may lead to widespread failures. Such adverse response was previously documented during several seismic events, such as the earthquakes of Niigata, Japan (1964, [7–10]), Dagupan City, Philippines (1990, [11–14]), Chi-Chi, Kocaeli, Turkey (1999, [15]), and Taiwan (1999, [16]).
The paper compares ordinary qualitative approaches generally used in liquefaction vulnerability assessment with a three-dimensional, advanced, and nonlinear finite element model results. The traditional methods simply consider the probability of a soil to be subjected to liquefaction without any quantitative considerations, and thus they are effective only in preliminary approaches. When more accurate predictions are required, nonlinear theories for counting the liquefaction behaviour must be adopted. In this paper, seismic effects during Tohoku earthquake are reproduced adopting a platform built up open-source computational platform OpenSees.
In the following sections, the state of the art about the most credited qualitative approaches, the employed computational formulation, and the calibration are described. The computational platform is then presented focusing on the adopted materials, boundary conditions, and analysis assumptions. Results of the conducted analysis are finally presented taking into account post-earthquake evaluations and eventually pre-earthquake considerations.
2. Urayasu Site
This study considers Urayasu site as a representative of Tohoku earthquake scenarios because of several aspects. First of all, even if the peak ground acceleration (PGA) value is ordinary, the registered liquefaction consequences (NIED [1], Towhata et al. [2], and Tokimatsu and Katsumata [6]) were extremely wide in terms of lateral spreading and settlements. Since these damages concentrated on residential ordinary buildings more than big cities, they were less affected by the influence of big and heavy structures. Consequently, free field conditions can be more representative than in urban areas. Finally, geological and geotechnical on-line database that referred to Urayasu site were detailed and well descripted.
Recordings and geotechnical data were deduced from National Research Institute for Earth Science and Disaster Prevention (NIED) [1] at Kyoshin Network K-Net. In particular, Figure 1 shows Urayasu (CHB008) Station (35.6537°N, 139.9023°E, and 379 km epicentre distance) recordings in terms of the east-west (EW) component acceleration time histories.
Geotechnical data taken as reference in this study such as values, shear velocity, and density versus depth diagrams and general characteristics layers are represented in Figure 2. The other parameters (i.e., friction angle and relative density) were derived by a calibration analysis shown in the next section.
3. Calibration Analysis
Urayasu stratification consists of two superficial sand and sand sandy layers (whose base is at about 6 m depth) upon a silt Pleistocene deposit. From an online database, it was possible to obtain the most representative geotechnical parameters.
Friction angles and relative densities were estimated from NSPT and effective vertical tension values through Schmertmann [17, 18] and Bazaraa [19] and Gibbs and Holtz [20] graphical and analytical correlations. The sand layers were characterized by and , taking into consideration that for the traditional nomenclature a sand layer contains more than 50% of sand material. Table 1 shows the estimated parameters. The silt deposit was modelled as a cohesive stratum characterized by medium density and cohesion parameters (Table 2).
Finally, water level is a significant parameter to take into account in liquefaction modelling, but it is also very difficult to determinate, because of its variable nature. In order to determine the most conservative assessment, water level was considered at the surface for all the models presented in this work.
Further detailed geotechnical definition will be needed to obtain more realistic quantitative response.
4. State of the Art
Traditional methods consist of defining liquefaction risk adopting very few (at least one) parameters that are able to represent complex liquefaction phenomenon. Generally, these methods follow two main directions: data from observed events and laboratory tests. Nowadays, a great number of existing criteria can be adopted: Sherif and Ishibashi [21], Youd and Perkins [22], and some others based on cyclic resistance ratio (CRR), such as Seed and Idriss [23, 24], Yegian and Whitman [25], Tokimatsu and Yoshimi [26], Iwasaki et al. [27, 28], Seed et al. [29], and Andrus and Stokoe [30]. More details are in [31, 32].
These empirical methods are generally built up with some significant parameters taken from the geological and geotechnical recognition, such as layer age, origin, and water depth. Some of these methods were transferred inside international codes such as Eurocode 8 [33], where liquefaction is verified considering various safety parameters adopting a simplified method based on empirical correlations between site measurements, data registered in previous earthquakes, and laboratory cycling test under controlled conditions. The next section details the application of these traditional methods to Urayasu City site.
5. Traditional Methods Response
Traditional methodologies verify the risk for the soil to be subjected to liquefaction without a direct control on excess pore pressure. They only refer to few descriptive parameters based on historic knowledge and on behaviour analysis results of laboratory cycling test under controlled condition. They are used only in preliminary studies allowing general considerations on liquefaction effects and deformation in qualitative terms. These methods instead are not sufficient if more accurate predictions are required, such as assessment of the displacements connected to lateral spreading, settlement predictions, or finally countermeasures design. In this regards, Yegian and Withman [25], Seed et al. [29] and Youd and Idriss [34] methods were applied to Urayasu site adopting the previous model definitions (Tables 1 and 2).
In particular, Yegian and Withman method is based on ILP index, and when this value is more than unity, as shown in Figure 3, the soil is potentially affected to liquefaction. The other methods define instead a safety coefficient () that has to be more than unity to preserve from liquefaction. Figure 3 shows that these methods verify liquefaction.
(a)
(b)
(c)
(d)
6. Computational Formulation
This section describes the computational approach. A method based on derivation of appropriate loading-unloading flow rules has been adopted in order to reproduce numerically the observed strong dilation tendency and the resulted increase in cyclic shear stiffness and strength. In particular, simulations were conducted using the open-source computational platform OpenSees [35]. This platform allows the implementation of a framework for saturated soil response as a two-phase material following the - (: soil skeleton displacement and : pore pressure) formulation of Chan [36] Zienkiewicz et al. [37]. Although based on these sophisticated theories, the model has the main advantage to be built up mainly with the most common used geotechnical parameters. Other parameters connected with the liquefaction mechanism can be obtained by assigned values calibrated on a big variety of realistic cases. In the study, two different models for cohesionless and cohesive behaviours were considered.
The first soil model is developed within the framework of multiyield surface plasticity theory for frictional cohesionless soils proposed firstly by Prevost [38]. In this model, emphasis is placed on controlling the magnitude of cycle-by-cycle permanent shear strain accumulation in clean medium to dense sands (Parra [39], Yang [40], and Yang and Elgamal [41]).
Cohesive soils are modelled as nonlinear hysteretic materials (Parra [39], Yang [40], and Yang and Elgamal [41]) with a Von Mises multisurface kinematic plasticity model. The focus is on reproduction of soil hysteretic elasto-plastic shear response (including permanent deformation). Constitutive model simulates monotonic or cyclic response of material whose shear behaviour is insensitive to confinement changes. Plasticity is exhibited only in the deviatoric stress-strain response and formulated based on the multisurface (nested surfaces) concept, with an associative flow rule. The volumetric stress-strain response is linear elastic and independent of the deviatoric response. Nonlinear shear stress-strain back bone curve is represented by one hyperbolic relation, defined by two material constants: low-strain shear modulus and ultimate shear strength.
The adopted parameters (especially for the cohesionless layers) were calibrated through an identification analysis in order to take into account the nonlinearity of liquefaction-induced behaviors. For more details, see Yang et al. [42], Elgamal et al. [43–45], and Forcellini and Tarantino [46].
7. Computational Framework
Simulations conducted in this paper apply the open-source computational platform OpenSees (Mazzoni et al. [35]) that can simulate performances of structural and geotechnical systems subjected to static and seismic loadings. In particular, this study adopts the interface OpenSeesPL (http://www.soilquake.net/openseespl/) originally calibrated for pile analyses and here modified in order to allow free field response studies. The interface simplifies the 3D spatial soil domain, boundary conditions, and input seismic excitation definition and also convenient postprocessing and graphical visualization of the analysis results including the deformed mesh and ground response time histories. OpenSees ability to simulate the real wave propagation adopting realistic boundaries is of particular importance and significance in order to realistically reproduce the damage scenarios.
The 3D m, 18 m high model built up with 2256 nodes and 1848 Brick-up hexahedral linear isoparametric 8-node elements is shown in Figure 4. Each node has 4 degrees of freedom (DOF), three DOFs for solid displacement (), and one DOF for fluid pressure () (Figure 5). This element is implemented for simulating dynamic response of solid-fluid fully coupled material, based on Biot’s theory of porous medium (for more details, see Lu et al. [47]).
The 3D model is assumed to have periodic boundary conditions, as to represent a small section of a presumably infinite (or at least very large) soil domain by allowing the energy imparted by the seismic event to be removed. First of all, this is obtained with the nodes on the base of the model free to displace in horizontal and fixed directions against the vertical translations. Secondly, periodic horizontal boundaries were chosen in order to ensure that free-field conditions really exist. In this regard, each node on the front boundary moves in the same way as the analogous node on the back boundary (horizontal and vertical directions).
In order to control the various steps of the problems and to manage the wise quantities of results, analysis was split into two consequent steps (gravity application and dynamic analysis itself). The first step ensures that the pore pressure distributions and effective stresses are appropriate for the site conditions prior to the application of a ground motion. Separate recorders were set up as to distinguish the gravity analysis from any other postgravity results. Nodal displacements, accelerations, and pore pressures are recorded along with the elemental stresses and strains at each of the nine Gauss points. In order to achieve hydrostatic pore pressure conditions, gravity application analysis divided itself into two parts where soil elements are considered to be firstly linear elastic and then plastic. The elastic portion of the gravity analysis is run as a transient analysis with very large time steps, thus simulating a static analysis. Therefore, the plastic portion of the gravity analysis is run using smaller time steps to aid in convergence. Finally, dynamic excitation analysis is developed using a transient analysis. For more details, see Elgamal et al. [45] and Forcellini and Tarantino [46].
8. Materials
A calibration procedure was fundamental in order to set up a three-layer (3 m + 3 m + 12 m) soil model, representative of 3D realistic strata. The superficial layers were modelled with a 3 m superficial sandy soil and a 3 m sand layer, adopting pressure dependent multiyield materials. Table 3 shows shear () and bulk modulus () at reference pressure (80 kPa), PT angle (), permeability (), contraction (), dilation ( and ), and liquefaction (, , and ) parameters for the superficial layers.
The silt deep layer was modelled as medium clay 12 m layer adopting the pressure independent multiyield model. Table 4 shows shear () and bulk modulus () at reference pressure (80 kPa), cohesion (), and permeability () parameters for the deep layer. For more details on the material definition, see Yang et al. [42], Elgamal et al. [44], and Lu et al. [47].
9. Numerical Response
This section shows numerical results in terms of excess pore pressure, longitudinal displacements, settlements time histories, and entire mesh deformation. These results confirm generally that liquefaction mainly affected the superficial layers.
In particular, Figure 6 shows the excess pore pressure time histories at 3.00 m, 6.00 m, and 12.00 m depths. Very low pore pressures resulted inside the base silt layer, while in the two superficial sand layers, values around 45 kPa for 6.00 m and 25 kPa for 3.00 m are reached and then dissipated after PGA value is reached. The result shows how the presence of these two layers (between 6.00 m and the surface) generates high pore pressure and thus liquefaction concentrated in these superficial layers. These results can be strengthened comparing 6.00 m, 3.00 m, and 0.00 m depth displacement time histories (Figure 7) that show how lateral spreading is totally due to the two superficial sand layers (between 6.00 m and the surface) that slide on the deepest one. In particular, lateral displacement at the top of base silt layer is zero, while maximum displacements at the top of the silty sand and at surface are around 12.20 cm and 16.40 cm, respectively. The entity of these displacements confirms that the site is strongly subjected to soil liquefaction and consequently to lateral spreading and is comparable to that observed in free field conditions.
(a)
(b)
(c)
(a)
(b)
(c)
Figure 8 shows the settlements time histories at the surface and at 3.00 m depth. These values (24.70 mm for the surface and 3.30 mm at 3 m depth) are representative of what is observed in Urayasu City just after the shake under free-field conditions.
(a)
(b)
Effects of liquefaction to the whole mesh are finally represented in Figures 9 and 10, where 3D and plan views of the entire mesh deformation in correspondence with PGA value are shown.
10. Conclusions
The paper compares ordinary qualitative approaches generally used in liquefaction vulnerability assessment with a three-dimensional, advanced, and nonlinear finite element model. Traditional methods simply verify liquefaction phenomenon existence without any quantitative response. Numerical simulation based on OpenSeesPL platform not only allows to define the susceptibly as all the traditional methods do, but it gives the possibility to predict numerical responses in order to evaluate the entities of damages.
In this regard, the obtained values are shown to be comparable to those observed in free field conditions, thus confirming Urayasu City site liquefaction vulnerability in terms of pore pressure built up, lateral spreading, and settlements.
The conducted analysis demonstrates the OpenSees high potentialities in performing appropriate numerical simulations for predicting liquefaction-induced lateral deformation. The results confirm the assumptions concerning the reproduction of the model in terms of calibration and numerical choices, such as boundary conditions, material definition and analysis choices.
In this regard, the response can quantify the performance and risk of liquefaction using metrics that are of immediate use for both pre-earthquake and post-earthquake risk assessment analyses. This kind of response can become very powerful if applied to soil-structure interaction studies. This will be object of further work.