Abstract
According to the investigation and restoration report by the Ministry of Land, Infrastructure, Transport and Tourism (MLIT, 2011), 1195 river embankments were damaged in Tohoku region during the 2011 Tohoku earthquake. The failures of the river embankments were typically due to the soil liquefaction of embankment fill. In the present study, a severely damaged river embankment along Naruse River was simulated by using a three-phase coupled finite element program, namely, COMVI2D-DY, which is developed to analyze large deformation behavior of partially saturated soils. In addition, the reconsolidation process after the earthquake was simulated. To reproduce the reconsolidation behavior, a cyclic elastoplastic constitutive model based on nonlinear kinematical hardening rule was modified by considering stiffness recovery during reconsolidation. From the analysis results, it could be concluded that the numerical method is able to reproduce the key characteristics of the actual damaged pattern; the embankment is heavily damaged and deformed largely towards the land side, and the settlement at the top of the embankment is 2.5 m. In addition, realistic simulation results can be obtained from the reconsolidation analysis.
1. Introduction
A river embankment is a large-scale earth structure that mainly functions as a measure for flood prevention. Though the geological stratigraphy of the river embankment varies from location to location, a particular common feature is that the ground water table within the embankment is very shallow, hence keeping the earth embankment close to full saturation. Such feature will make the embankments susceptible to soil liquefaction, especially for those which are filled with granular soils.
In Japan, the seismic performance of the river embankment has been inspected in the last decades [1] since the significant subsidence of the Yodo River levee in Osaka was reported by the 1995 Kobe earthquake. The vulnerability of numerous river embankments, however, remains unknown until the occurrence of embankment failure due to soil liquefaction induced by large earthquakes. For instance, more than 2000 river embankments were damaged during the 2011 off-the-Pacific-Coast Tohoku earthquake. According to the investigation and restoration report by the Ministry of Land, Infrastructure, Transport and Tourism [2], 1195 river embankments were damaged in the Tohoku region with the total length of approximately 100,000 meters. Liquefaction of river embankments and the failure mechanism during the 2011 Tohoku earthquake were discussed and summarized by some researchers (e.g., [1, 3–5]). The failures of the river embankments were mainly due to the soil liquefaction of foundation soil or embankment fill. Figures 1 and 2 show heavily damaged embankments during the earthquake.
In order to prevent future losses, it is necessary to understand the failure mechanism of river embankments due to soil liquefaction. On the other hand, the restoration time for damaged river embankments is usually limited, as they are to be reinstated before the raining reason. It is hence common that the embankment was repaired without the knowledge of its stress status or hydraulic characteristics. This would result in potential risk of damage in the future. Hence, a reconsolidation mechanism after liquefaction is as important as a soil liquefaction mechanism for a river embankment.
Over the past few decades, multiphase coupled numerical analysis methods have been proposed for partially saturated soils (e.g., [6, 7]). Numerical simulations of embankment liquefaction have been conducted using various modelling methods in conjunction with a variety of constitutive models. Some (e.g., [8–10]) adopted infinitesimal strain theory to simulate seismic response of embankments. The infinitesimal strain theory is capable of capturing some important characteristics of soil liquefaction at the element level, such as the increase in pore water pressure and the loss of effective stress, but can hardly replicate the embankment deformation pattern, especially under the extreme ground motion. In addition, the behavior of partially saturated soil is not considered in most, and the embankment fill is generally treated as a dry material. To analyze the soil behavior within the embankment, however, the descriptions of partially saturated soil, such as the soil-water characteristics and the constitutive model for unsaturated soil, need to be introduced. Recently, Yoshikawa et al. [11] conducted a numerical analysis on the seismic behavior of an unsaturated embankment focusing on the effect of the groundwater level. There are, however, few numerical examples to reproduce the actual damage of river embankment.
In the present study, the dynamic behavior of a severely damaged river embankment along Naruse River during the 2011 Tohoku earthquake was simulated by using a three-phase coupled finite element program, namely, COMVI2D-DY [12, 13]. The program COMVI2D-DY was developed to analyze large deformation behavior of partially saturated soils under dynamic loading. The simulations are compared with the observed deformation pattern of the river embankment. Failure mechanisms due to liquefaction of embankment fill including but not limited to the pore water pressure distribution and the effective stress distribution can be clarified. In addition, reconsolidation analysis was carried out as a continuation of the dynamic one to understand the mechanical as well as the hydraulic characteristics of the damaged embankment before restoration. To reproduce the reconsolidation behavior, a cyclic elastoplastic constitutive model based on the nonlinear kinematical hardening rule [14] was modified by considering stiffness recovery during reconsolidation, and the behavior is discussed.
2. Governing Equations for Large Deformation Dynamic Analysis of Multiphase Geomaterials
A finite element program, namely, COMVI2D-DY [13], was developed to simulate large deformation of partially saturated soil under dynamic loading. The governing equations of dynamic analysis for multiphase geomaterials, that is, air-water-soil three-phase materials, are described within the framework of theory of porous media. In the theory, the governing equations are described for the multiphase mixture, that is, an immiscible mixture of solid and fluids [15]. The u-p formulation based on the updated Lagrangian method has been employed along with the Jaumann rate of Cauchy stress tensor for the weak form of the equilibrium equation.
An elastoplastic constitutive model [14, 16] is used to reproduce cyclic plastic behavior and liquefaction behavior, that is, the reduction of the mean effective stress of the saturated sandy soil. In the model, a nonassociativity parameter that controls the nonassociative behavior of sandy soils is incorporated. The formulation of the elastoplastic constitutive model is given in the following section. For the partially saturated sand, an extended elastoplastic model considering the effect of suction [17] is used, in which the so-called skeleton stress which provides a natural application of the two-phase mixture theory to the three-phase mixture (e.g., [18–21]) is the basic stress variable. In addition, the suction is included in the overconsolidation boundary surface and the yield surface to describe the bonding effect. An elastoviscoplastic constitutive model by Kimoto et al. [22] was used for clayey soils which can reproduce the dynamic behavior of clayey soils which exhibit plastic behavior as well as the rate-dependent behaviors. A van Genuchten-type model [23] is used as a soil-water characteristic curve.
The weak form of the continuity equation for water and gas, the conservation of momentum for the whole mixture is discretized in space by the finite element method. Assuming two-dimensional plane strain condition, an 8-node isoparametric element is used for the displacement, velocity, and acceleration of the solid skeleton, and a 4-node isoparametric element is used for pore water pressure and pore gas pressure. For the time discretization, Newmark’s β method is adopted. The detailed descriptions for the dynamic finite element formulations are shown by Shahbodagh Khan [13]. Note that the pore-air pressure is assumed to be zero in the following simulations since the effect of air pressure was found to be negligible under drained conditions at the surface boundary.
3. Elastoplastic Constitutive Model
3.1. Formulation of the Model
Oka and Kimoto [14, 16] have developed a cyclic elastoplastic constitutive model based on the kinematical hardening theory by extending the previous model by Oka et al. [24]. In the new model, the yield function includes the changes in the stress ratio and the mean effective stress. A nonassociativity parameter that controls the nonassociative behavior of sandy soils, which is important for the simulation of unstable behavior such as liquefaction, is incorporated. In this section, the formulation of the model is briefly shown. The detailed explanations for the model are shown in Oka and Kimoto [14, 16].
The overconsolidation boundary which defines the boundary between the normally consolidated region and the overconsolidated region is given in
In this equation, indicates a normally consolidated region, whereas indicates an overconsolidated region. Relative stress ratio is given in where is the mean effective stress and denotes the stress ratio tensor. The hardening parameter controls the size of the overconsolidation boundary surface. is the value of the stress ratio at the transformation phase, at which point, the soil transforms from compression to dilation. The initial value of is equal to the stress ratio tensor after the end of consolidation, and follows and evolutional equation as where is the accumulated plastic shear strain from the initial state, is the plastic deviatoric strain increment, and is a material constant which controls the disappearance of anisotropy.
The static yield surface is given by a function of the nonlinear kinematic hardening rule of the stress ratio, mean effective stress, and plastic volumetric strain, in where is in the unit of stress which is kilopascal (kPa). is the kinematic hardening parameter that associates with the change of plastic volumetric strain. Parameter is the nonassociativity material parameter which is ranging from zero to one; for associated flow rule and for nonassociated flow rule. It should be noted that is assumed before the stress ratio reaches the phase transformation line and during consolidation after liquefaction. is the back stress parameters with the same dimension as the stress ratio . The nonlinear hardening rule, as an evolutional equation of , is given by where and are material parameters related to the stress ratio at failure, , and the nondimensional initial plastic shear modulus, , as follows:
In addition, follows an equation as where is the lower boundary of . Parameter controls the rate of hardening, and is the maximum value of , which is expressed in where is a hardening parameter before the stress ratio reaches the phase transformation stress ratio , is the reference plastic strain, and is accumulative plastic shear strain after phase transformation. Note that the material parameter is set to zero in this research.
The kinematical hardening associated with the change in the mean effective stress is given by where , , and are material constants and is the volumetric plastic increment.
The plastic potential function is given by where is a constant and is the dilatancy parameter whose definition is distinguished between the normally consolidated region and the overconsolidated region, as given in where is a material parameter and is the mean effective stress at the intersection between the surface which passed through the current stress and has the same shape as , and the anisotropic stress axis. The overconsolidated boundary surface, yield function, and plastic potential function for an isotropic overconsolidated state are illustrated in Figure 3.
Oka and Kimoto [14] used Equation (12) to account for the stress and strain dependency of the elastic shear modulus: where , , and are material parameters and the initial elastic shear modulus. denotes the accumulative plastic deviatoric strain after the stress ratio reaches the phase transformation , and is the reference strain which can be obtained by fitting to liquefaction resistance curve.
3.2. Parametric Study on
The effect of nonassociativity parameter in the cyclic elastoplastic model under undrained simple shear conditions was presented in Oka and Kimoto [16] by simulating the behavior of loose Toyoura silica sand. They showed that the parameter controls the minimum mean effective stress during cyclic loading, namely, the minimum mean effective stress is smaller for the case with smaller . In this section, a parametric study is performed by using COMVI2D-DY with a 1-D horizontally layered model to demonstrate the effect of on the dynamic behavior and the reconsolidation behavior. The 1-D model consists of 10 elements with a total height of 30 m and a width of 3 m as shown in Figure 4. The bottom boundary of the model is perfectly fixed, and the equal nodal displacements on both sides are assumed at the same depth. A pervious boundary with a zero static pore water pressure is set at the top of the model while an impervious boundary is set for the other boundaries. The ground motion of the 2011 Tohoku earthquake recorded at a depth of 80 m at the Tajiri earthquake observation station in Miyagi Prefecture [25] is used. The acceleration profile is shown in Figure 5. The dynamic loading lapses after 200 seconds; the modelled soil column is then allowed to consolidate for 10 days. Material parameters used in this parametric are shown in Table 1 in the right column, which is determined for Naruse River sand in Miyagi Prefecture.
The effect of nonassociativity parameter is revealed by comparing analysis Case 1: , Case 2: , and Case 3: . Figures 6(a) and 6(b) show the time history of the pore water pressure and the mean effective stress of element 8. It can be observed that has almost no effect on the excess pore water pressure generation during dynamic loading. However, a higher value of leads to faster dissipation of the excess pore water pressure during reconsolidation, as well as a faster recovery of the mean effective stress. The final surface settlements of case 1, case 2, and case 3 are 0.68 m, 1.11 m, and 1.18 m, respectively, as shown in Figure 6(c). It appears that the value of has significant influence on consolidation settlement. It is worth noting that associate flow rule () is assumed during consolidation analysis for all the three cases. Therefore, the difference in consolidation settlement between each case is due to the combination effect of the value and accumulative plastic strain during dynamic analysis.
(a) Effect of on time history of pore water pressure of element 8
(b) Effect of on time history of mean effective stress of element 8
(c) Effect of on time history of surface settlement
(d) Stress path of element 8 for Naruse sand with Case 1 and Case 3
Moreover, controls the reduction of mean effective stress when soil liquefies, together with the nonlinear kinematic hardening parameter . Therefore, the minimum mean effective stress is larger in Case 3 than in Case 1 as shown in Figure 6(d). This tendency agrees with the results shown in Oka and Kimoto [16]. Besides, the response of stress due to ground acceleration is slightly smaller for larger . Such a mechanism is important for not only the realistic soil behavior simulation but also the stability of numerical calculation. Furthermore, larger increase in stress ratio in Case 3 is found during consolidation than that in Case 1. This is due to the higher stiffness of soil in Case 3 which would result in a higher response stress.
From the parametric study, it was found that the parameter has a significant influence on consolidation settlement, that is, a larger value of gives smaller settlement for this set of parameters. It should be noted, however, that the model with strong nonassociativity with a small value of has a potential for simulating liquefaction while the model with the associated flow rule does not as mentioned in Oka and Kimoto [16]. In other words, a large value of may not be appropriate to simulate liquefaction behavior.
3.3. Modification of the Model considering Stiffness Recovery
In the elastoplastic constitutive model shown in the previous section, the reduction of the elastic stiffness and the plastic stiffness with increasing deviatoric strain was given in (8) and (12). In the following, the elastic and plastic stiffness recovery during reconsolidation is introduced in the model.
The elastic shear modulus is extended to incorporate the hardening effect of volumetric strain compression during reconsolidation in the form of in which is a material parameter that represents the value of final elastic shear modulus at the end of consolidation and is a material parameter to control the recovery rate of elastic shear modulus during consolidation. The accumulative volumetric compression after dynamic loading, namely, during reconsolidation, is denoted by . Equation (13) separates the hardening component attained from volumetric compression, that is, consolidation, from that attained during dynamic shearing. In addition, it is assumed that is not affected by the change of plastic deviatoric strain during consolidation. In other words, is constant.
Similar to the elastic shear modulus, the hardening term is considered in the plastic shear modulus during consolidation which is associated with the volumetric compression: where is the material parameter and equals the value of recovery at the end of consolidation. Ye et al. [26] experimentally investigated the stiffness recovery of liquefied sandy ground and showed that the stiffness after earthquake motion was the same or larger than that before the earthquake. In the present study, and are assumed.
A comparison between the original elastoplastic model [24] in which the yield function is given by the stress ratio, and the extended elastoplastic model with stiffness recovery will be presented in the following. The exercise is mainly to demonstrate the stiffness recovery ability of the modified model with modifications in the current research. The same model configuration, drainage conditions, boundary conditions, and loading conditions as those used in the previous section are adopted. Material parameters are shown in Table 1. is chosen for the case of strong nonassociativity, similar to the cyclic elastoplastic model. The hardening parameter is set to be . It was widely understood that liquefied soil would experience a low-rigidity period before the stiffness recovery. Therefore, the onset of consolidation in the case of using the extended model is set to be 4000 seconds, which correspond to the time lapse before the excess pore water pressure starts to dissipate in the case of the original model.
Simulation results are shown in Figures 7(a)–7(g) for element 8. It is observed that the cyclic elastoplastic model [24] experienced higher excess pore water pressure during dynamic loading than the extended model. This is mainly due to the combination effect of and the kinematic hardening term. Figure 7(b) shows the comparison of the mean effective stress of element 8, in which the elastoplastic model experienced lower mean effective stress. The effect of volumetric compression hardening on stiffness recovery is depicted in Figures 7(c) and 7(d). Figure 7(c) shows the recovery of elastic shear modulus, whereas the recovery of the kinematic hardening parameter , which is related to soil plastic stiffness, is shown in Figure 7(d). The stiffness recovers effectively and ends with a higher final stiffness than the initial one. Note that the amount of stiffness recovery at the end of consolidation is an input parameter and should be determined from experimental data. However, due to the lack of information, they are assumed to be equal to the initial stiffness value. Figures 7(e) and 7(f) show the effective stress path of element 8 for the original elastoplastic model and extended elastoplastic model, respectively. The minimum mean effective stress is smaller in the original model. Significant distinction is observed in Figure 7(g)—the total settlement of the case with the elastoplastic model is 0.67 m whereas that of the case with the extended model is 0.29 m. This is due to the effect of kinematic hardening and stiffness during consolidation. Such a result is reasonable because it is not realistic for sandy soil to have much larger settlement due to consolidation than during liquefaction.
(a) Comparison of time history of pore water pressure
(b) Comparison of time history of mean effective stress
(c) Comparison of time history of elastic shear modulus
(d) Comparison of time history of kinematic hardening parameter
(e) Stress path of element 8 for cyclic elastoplastic model
(f) Stress path of element 8 for extended cyclic elastoplastic model
(g) Comparison of time history of settlement
4. Numerical Analysis of River Embankment during Earthquake and Reconsolidation
4.1. Model Description
The target river embankment is located at Shimonakanome upper stream of Naruse River, 30.3 km from the river mouth. The location of the embankment is shown in Figure 8. This embankment was found heavily damaged after the 2011 Tohoku earthquake. According to the investigation report [27], boiling sand was observed in the longitudinal cracks as shown in Figure 9. In addition, the presence of ground water in the embankment was confirmed. Therefore, it was concluded by MLIT [2] that the failure of river embankment is caused by liquefaction of the embankment fill.
A cross section of the river embankment including its geometry, geological stratigraphy, and the location of the water table is shown in Figure 10. The liquefied soil is identified as clayey sand (B1s) with a maximum thickness of 2.9 m. According to the geological map of Naruse River, this river embankment was a natural levee [28] which was indicated as “old embankment” in Figure 10. The old embankment was formed by sandy soil (B1c and B1s) and then enlarged twice using sand and clay material. The location of ground water table is indicated by the blue line, which was measured on site one month after the earthquake. The deformation shape is illustrated in Figure 10 as well. It could be observed that large sliding occurred towards the land side, after liquefaction.
Geometry of the finite element model, the mesh, and the boundary conditions are depicted in Figure 11. The numerical simulation was conducted under 2-D plane strain condition. A full-scale model was considered with an engineering datum of 50 m below ground level. To reduce the side boundary effects, semi-infinite elements with equidisplacement conditions were used on the left and right ends, each being 100 m long. Drainage boundary is set at ground level where the initial pore water pressure is zero. A no-water-flow boundary is prescribed along the surface of the embankment to avoid unrealistic water inflow due to suction. Regarding the mesh size in the vertical direction, it is generally determined to be 1/5 of the minimum wavelength in dynamic FEM analysis. For example, when the shear wave velocity is 100 m/s and the target maximum frequency is 10 Hz, the minimum wavelength is 10 m; thus, the required mesh size is 2 m. In the present analysis, the vertical mesh size of the embankment is determined to be around 0.5 m considering the change of the shear wave velocity and the frequency after liquefaction.
Ground acceleration of the Tohoku earthquake recorded at a depth of 80 m at Tajiri station in Miyagi Prefecture shown in Figure 5 was used during dynamic analysis, which is the same as that was used in the 1-D dynamic-reconsolidation analysis. The rigid base is used so that the bottom boundary of the model is fixed. The Tajiri earthquake motion observation station is located approximately 8.7 km away from the river embankment at Shimonakanome upper stream of Naruse River.
4.2. Material Parameters
In this analysis, sandy soil is modeled by the extended elastoplastic model whereas the elastoviscoplastic model is used for clayey soil. Soil physical tests were carried out and presented in the investigation report of the damaged river embankment in the Tohoku region [2]. It was found that fine content of the embankment fill is generally close to 50%; hence, the extended elastoplastic model is adopted for all the embankment fills for the sake of numerical simplicity. For the partially saturated region above the water table, namely, the B1c, B2c, and B3s layers, an extended elastoplastic model considering the effect of suction is used, in which the skeleton stress is the basic stress variable and the suction is included in the overconsolidation boundary surface and the yield surface to describe the bonding effect. The embankment is situated on a thick layer of soft clay, which is simulated by the elastoviscoplastic model. Note that the foundation soil is assumed to be homogenous in this analysis.
The cyclic behavior of soil can be represented by a liquefaction resistance curve. Two cyclic triaxial tests were conducted for this river embankment [25]. Location of the soil samples are shown in Figure 10. Due to inadequate experimental results, cyclic behavior of soil B2c and B1c is considered as the same as that of B3s and B1s. Soil test simulation was carried out, and the material parameters are determined to reproduce the results from cyclic triaxial test. The cyclic triaxial test for B3s and B1s and their numerical behavior under cyclic load are plotted in Figure 12. Since the cyclic triaxial test was not carried out for the clayey foundation, material parameters for Torishima clay [13] are adopted since it is a soft alluvium deposit as well. The liquefaction resistance curve of Torishima clay is shown in Figure 13. Material parameters are shown in Table 2, and the suction parameters for the partially saturated region are shown in Table 3.
In COMVI2D-DY, the degree of water saturation depends on the soil-water characteristic curve (SWCC), adopting van Genuchten’s [23] model. In Carsel and Parrish [29], typical values of SWCC parameters, and , in van Genuchten’s model were proposed according to soil descriptions and their permeability. In this simulation, and , for silt loam, were chosen by considering the soil type and average permeability of the embankment fill. Hydraulic properties and material parameters for the effect of suction for unsaturated soil are listed in Table 3. Since the water table is located above the foundation soil, it is not necessary to consider the soil-water characteristic for soft clay. The SWCC used in the dynamic analysis is shown in Figure 14. The initial degree of saturation above the water table is assumed to be 0.6, which corresponds to an initial suction of 14.5 kPa.
In this analysis, an initial stiffness-dependent type of Rayleigh damping is adopted with an attenuation constant of 0.01. Time is discretized by Newmark’s method where and are 0.3025 and 0.6, respectively. The time increment for the dynamic analysis was set to be 0.002 second.
Before conducting the dynamic analysis, the initial effective stress was generated by a static finite element analysis with an elastoperfect plastic model of which the failure surface follows the Drucker-Prager model.
4.3. Dynamic Behavior
The comparison between the original shape of the river embankment, its observed deformed shape after the earthquake, and the deformed mesh from the dynamic analysis at the end of ground acceleration, 200 sec, is shown in Figure 15. It could be observed that the embankment is heavily damaged at the end of ground shaking. The key characteristics of the damaged pattern can be reproduced by dynamic analysis. First, the embankment deformed largely towards the land side. Second, part of the embankment on the river side remains undamaged. It appears that the actual deformation on the land side is much larger than that from the simulation. This large lateral spreading of the river embankment is not only due to the collapse of earth structure but also due to the sliding of soil. Such phenomenon is difficult to simulate by finite element analysis since this would involve disconnection between elements or soil particles. In addition, the amount of the simulated deformed soil cannot be compared to that which was observed on-site, since the ground acceleration used in this analysis is not at the exact location of the river embankment.
Time history of the horizontal and the vertical displacements at node 1201 and node 2193 is shown in Figures 16 and 17, respectively. The locations of node 1201 and node 2193 are depicted in Figure 18. It could be observed that the maximum lateral movements towards the land side and settlement at the top of the embankment are 5.02 m and 2.50 m, respectively. The simulated settlement at the top of the embankment is close to the observed settlement after the earthquake which is approximately 3 m. However, a lateral displacement of 12 m was observed on site which is much larger than the simulated one. This result agrees with the previous conclusion that the sliding of earth is difficult to simulate by finite element method. In addition, it can be seen that both large horizontal and vertical displacements occurred from 30 sec to 150 sec which corresponds to the time period from the start of the first main shock till the end of the second main shock.
Figures 19 and 20 show the acceleration response for node 1201 and node 2193, respectively. Both nodes exhibit acceleration amplification. The peak ground acceleration of the Tohoku earthquake occurs at 33.9 second and 81.2 second. An average delay of 2 seconds of the peak response can be observed in both nodes. Moreover, although the shaking reduces significantly after 150 seconds, the amplification of acceleration for both nodes remains large. This might be contributed by the viscosity of the soft foundation clay.
Failure mechanism of the river embankment can be clarified by examining the change of pore water pressure, change of stress, change of strain, and so on The increase of the pore water pressure results in reduction of the effective stress, which is a direct indication of soil liquefaction. Time history of the pore water pressure and the mean effective stress of element 681 and element 865 are plotted in Figures 21 and 22, respectively. Element 681 and element 865 are located at the center of the liquefied soil layer, as shown in Figure 18. It can be seen that the excess pore water pressure generates rapidly from about 30 seconds, the first main shock, and reaches maximum even before the second main shock. The mean effective stress experienced rapid decrease as well, corresponding to the time when a large amount of the excess pore water pressure is generated. Therefore, element 681 and element 865 are liquefied after the first main shock.
The distribution of the pore water pressure, the mean effective stress, and the skeleton stress decreasing ratio are shown in Figures 23, 24, and 25, respectively. It can be observed that the excess pore water pressure is generated within the embankment, which leads to the reduction of the mean effective stress to a minimum value of 1.9 kPa at 60 seconds. The mean skeleton stress decreasing ratio (ESDR) is defined by
According to Figure 25, the minimum value of ESDR in the embankment is 90% at 60 seconds, indicating a complete liquefaction after the first mean shock.
Distribution of the accumulative deviatoric plastic/viscoplastic strain (GAMP) is shown in Figure 26. Although GAMP appears to be small at 60 s, the maximum GAMP is in fact 52%. As the simulation continues, a clear shear band was developed along the base of the liquefied soil layer. At the end of the dynamic analysis, the maximum value of GAMP is 224% that causes large deformation of the embankment.
The distribution of the suction and distribution of the degree of saturation are shown in Figures 27 and 28, respectively. The initial maximum suction is set to be 14.5 kPa, as a result of assuming that the initial degree of saturation is 0.6. The maximum value of suction at the end of consolidation is 19.1 kPa. Therefore, in this analysis, the change of suction is not significant. It is observed that during the shaking, the unsaturated soil on the land side of the embankment increases whereas that on the river side decreases. This is mainly due to the deformation of the embankment. In addition, the increase in degree of saturation is also due to the shorter drainage length on the river side of the embankment that was caused by large settlement of the embankment.
4.4. Reconsolidation Behavior
Although the dynamic behavior and the reconsolidation behavior of liquefied river embankment are discussed separately, the reconsolidation analysis is a continuation of the dynamic simulation by setting ground acceleration to zero. Therefore, material parameters, boundary conditions, and numerical conditions remain unchanged from the dynamic analysis. Note that the reconsolidation analysis starts at 4000 seconds which is the same as the onset time in 1-D analysis, in order to retain the stiff low-rigidity zone of liquefied soil as well as to avoid any numerical instability. The time step for dynamic analysis was 0.002 second, which is too small for consolidation analysis. Hence, time steps are gradually increased during consolidation analysis. The time history of the time steps is shown in Table 4.
The results from reconsolidation analysis are shown in Figures 29, 30, 31, 32, 33, 34, 35, and 36. The same elements, element 681 and element 865, are chosen for illustrating the pore water pressure dissipation during reconsolidation. From Figure 29, it could be observed that the excess pore water pressure dissipates immediately after the beginning of reconsolidation, from 4000 second. Meanwhile, the mean effective stress is effective in recovering as excess pore water pressure dissipates, as shown in Figure 30. The excess pore water pressure is fully dissipated after approximately 400,000 seconds (110 hours). Fluctuation of pore water pressure and mean effective stress are found at about 60,000 seconds. One shall note that unlike one-dimensional analysis, the behavior of the elements is affected by global behavior of the river embankment and the interaction with other soil layers. The stress paths of element 681 and element 865 during dynamic and reconsolidation analysis are shown in Figure 31. It could be observed that the mean effective stresses for both elements are not completely reduced to zero after liquefaction. This is due to the combination effect by the nonassociativity parameter and the nonlinear kinematic hardening parameter .
(a) Element 681
(b) Element 865
Figures 32, 33, and 34 depict the change of distribution of pore water pressure, mean effective stress, and ESDR. The distribution of pore water pressure was smooth at 11 hours, but the excess pore water pressure was not fully dissipated. A full dissipation can be obtained after 110 hours. The distribution of ESDR also shows that the mean effective stress is effectively recovered at 110 hours. It was also found that the consolidation analysis has almost no effect on the distribution of suction and degree of water saturation.
As mentioned earlier, soil stiffness was described by the elastic shear modulus in (8). Rapid stiffness recovery occurred almost immediately after the onset of reconsolidation, as shown in Figure 35. Note that the stiffness would not recover until the occurrence of volumetric compression during reconsolidation. In addition, the softening effect due to the increase in accumulative deviatoric strain is seized during reconsolidation. The rate of stiffness recovery is mainly governed by the input parameter , which should be chosen depending on experimental results. In this analysis, is set to be 2000.
Vertical displacement during dynamic and reconsolidation analysis is shown in Figure 36. It can be seen that the vertical displacement rate decreased drastically after 200 seconds when the shaking was stopped. Consolidation settlement continues slowly for node 1201 and node 2193, at almost the same rate. The consolidation settlement was not only contributed to by the embankment fill but also by the foundation clay.
5. Conclusions
In the present study, first, an extended cyclic elastoplastic constitutive model [14] was introduced into a three-phase coupled finite element program COMVI2D-DY. This model is then modified by incorporating volumetric compression hardening during reconsolidation. The cyclic-reconsolidation behavior of the extended model is studied and examined by one-dimensional analysis. It can be concluded that stiffness recovery during consolidation is effectively reflected on the elastic and plastic shear stiffness, which results in significant reduction of the consolidation settlement. Second, a severely damaged river embankment along the Naruse River in Tohoku District during the 2011 off-the-Pacific-Coast Tohoku earthquake is simulated to validate the numerical method of COMVI2D-DY incorporating the behaviors of the cyclic elastoviscoplastic model and the extended cyclic elastoplastic model. The validation exercise is based upon comparing the monitored deformation pattern of the river embankment with that from simulation. From the simulation results, it could be concluded that the numerical method is able to reproduce the key characteristics of the actual damaged pattern.
Data Availability
The data used to support the findings of this study are available from the corresponding author upon request.
Conflicts of Interest
The authors declare that there is no conflict of interest regarding the publication of this paper.
Acknowledgments
The authors acknowledge the National Research Institute for Earth Science and Disaster Prevention (NIED) for providing KiK-net strong-motion data.