Nonlinear Dynamic Analysis Adopting Effective Stress Approach of an Embankment Involving Liquefaction Potential

. Stability of an embankment under earthquake loads is challenging in the process of analysis and design. Some embankment design consist of saturated granular material that is potential to liquefaction. Earthquake loads to the embankment under this conditions is one of major cause of embankment failure. Seismic performance involving stress-deformations and excess-pore-water pressure was evaluated in this paper. The evaluation adopts effective stress approach with non-linear elasto-plastic constitutive model. Numerical simulations through parametric studies were performed to estimate minimum density and embankment height efficiently to tolerate lateral displacements due to liquefaction. A number of parametric analyses were performed to investigate the relationships among relative densities of sand, ground accelerations, embankment height to excess-pore-pressure and lateral displacement of the embankment. The liquefaction analysis is conducted numerically using a finite difference method FLAC Dynamic 2D software adopting Finn-Byrne constitutive model.


Background
Past earthquakes confirmed that the major cause of severe damages to embankments is the liquefaction of the soil foundation.Flow liquefaction and cyclic mobility deserve special attention because of the large deformations that accompany these two phenomena.The performance of embankment subjected to seismic action can be evaluated through different approaches including force-based pseudo-static methods, displacement-based sliding block methods and dynamic analysis (Sherard, 1967;Seed et al., 1975).Gazetas (1987) outlined important features, advantages and limitations of the methods.Pseudo-static approach, the most common methods is widely used in engineering practice to assess the seismic stability of embankment.This approach is quite simplistic by which the complex aspects of seismic behaviour are represented in terms of static forces and the embankment stability is expressed in terms of an overall factor of safety.The response of embankment to earthquake may be related to several factors such as embankment geometry, mechanical properties of construction soil materials, distributions of pre-seismic stresses and pore-waterpressures, and input motion characteristics.Most of these factors are partially or totally neglected by the approaches traditionally adopted for assessing the seismic safety of embankment.For instance, some earthquake parameters such as frequency content and duration which significantly affect the soil response are neglected in the pseudo-static approach (Ambraseys, 1960).
The development of a better understanding of the fundamental mechanisms leading to the generation of pore pressures during undrained cyclic loading of sands (Martin et al., 1975), resulted in the development of an alternative analytical approach based on an effective stress model (Finn et al., 1977).In this approach, pore pressure increases are coupled to dynamic response solutions, enabling the complete time history of pore pressure increases to be computed during an earthquake.This method also allows the effects of soil stiffness degradation resulting from pore pressure increases to be reflected in the dynamic response solutions.Furthermore, the effects of pore pressure redistribution and dissipation can also be taken into account.
The seismic effective stress analysis is one of the most advanced numerical analyses used in geotechnical engineering.It allows to simulate accurately complex dynamic behaviour of soils during earthquakes including soil non-linearity, real earthquake records, rapid development of pore water pressures, soil liquefaction, and their effects on foundations/structures.
Recently, the seismic effective stress analysis has been adopted as a state-of-the-art/practice for evaluation of liquefaction problems including assessment of liquefaction-induced displacements, effects on structures and effectiveness of countermeasures against liquefaction (TFR,2007;ISO,2005) In this regard, parametric studies are conducted to identify the effects of embankment height, input motion characteristics, soil behaviour on the seismic response of embankment.It should be mentioned that this study does not consider the fluid-skeleton interaction, which may have significant effects on the seismic response of embankment.The liquefaction analysis is conducted numerically by using a finite difference method FLAC Dynamic 2D software adopting Finn-Byrne constitutive model.

Methodology of Effective Stress Analysis
The seismic effective stress analysis involves the following main steps: 1. Definition of the numerical model 2. Determination of parameters of the constitutive model

Dynamic analysis and interpretation of results
Firstly, the numerical model is defined by selecting appropriate element types, dimensions of the model, mesh size, boundary conditions and initial stress state of the soil.The geometry of the structure and stratification of the soil deposit needs to follow the basic rules of a good numerical analysis in order to define a numerical model.The problem selected here is the simplified representation of typical embankment consist sand material over clay layer.Furthermore, the selection of appropriate boundary conditions along end-boundaries or soil-foundationstructure interfaces is critically important in order to allow development of unconstrained response and deformation/failure modes.Similarly, an initial stress analysis is required to determine gravity-induced stresses in the soil and account for their effects on the stress-strain behaviour and liquefaction resistance of soils.
In the second step, parameters of the constitutive model need to be determined using data from field investigations and results from laboratory tests on soil samples.The types and number of laboratory tests required for the parameters of the constitutive model may vary significantly and are model-dependent.It is important to note that whereas most of the constitutive model parameters can be directly evaluated from data obtained from laboratory tests and in-situ investigations, there are always some parameters that are determined through a calibration process often referred to as element test simulations.Bes-fit values of the model parameters are identified through this process simulations of observation soil behaviour in laboratory tests.
In the final step, an acceleration time history is selected as an input ground motion, which is typically defined as a base excitation for the model.Considering the geometry of the problem and anticipated behaviour, numerical parameters such as computational time increment, integration scheme and numerical damping are also adopted, and the dynamic effective stress analysis is then executed.The analysis is quite demanding on the user in all steps including the final stages of postprocessing and interpretation of results since it requires an in-depth understanding of the phenomena considered, constitutive model used and particular features of the numerical procedures adopted in the analysis.Fig. 1 summarizes the main steps and procedures involved in a seismic effective stress analysis (Cubrinovski, 2011).

Constitutive Model
The ability of the seismic effective stress analysis to accurately simulate soil behaviour during earthquakes essentially depends on the ability of the constitutive model to represent real soil behaviour.In the case of soil liquefaction, this is a demanding task due to the complexity of the phenomena considered.Large and abrupt variation in soil stiffness, significant reduction in strength, and post-liquefaction soil behaviour are unpredictable parameters.
While recognizing, the art of modelling is really about selecting and focusing on the most important aspects of the problem considered.Thus a good constitutive model will be tailored to represent the key aspects of soil behaviour.
In the static analyses, the Mohr Coulomb model was used to represent the constitutive behaviour of all the material zones of the embankment.For the dynamic analyses, the Finn-Byrne model was chosen to simulate the constitutive behaviour of the granular material zones, and the Mohr Coulomb model was applied to base layer which are assumed to be non-liquefied.A detailed description of the two selected constitutive models is presented in the following sections.
The Mohr Coulomb (MC) model is a simple linear elastic perfectly-plastic model which requires five input parameters.These input parameters are Young's modulus, Poisson's ratio, friction angle, cohesion, and dilatancy angle.These parameters can be obtained from basic soil tests.In the finite difference program FLAC, some additional material parameters such as unit weight and hydraulic conductivities are also needed.To be remembered that for the dynamic analyses the MC model is just utilized for the non-liquefiable zones.
The value of Poisson's ratio is assumed to be 0.3 for all the materials.This value of Poisson's ratio is considered to be suitable for this type of analysis.
The value of the dilatancy angle is assumed to be zero for all material zones in this analysis.This is a convenient assumption, the unrealistically high negative pore pressures may occur due to the use of a positive value of the dilatancy angle, whereas unreasonably large positive pore pressures may develop if the value of the dilatancy angle is negative.(Brinkgreve, 2011) Fig. 1.Main steps and procedures in seismic effective stress analysis (Cubrinovski, 2011) Finn-Byrne is an effective stress elastic-plastic model which is capable of simulating the liquefaction behaviour of sands and silty sands under seismic loading.Using this model, it is possible to calculate pore water pressure generation by calculating irrecoverable volumetric strains during dynamic analysis.The void ratio in this model is supposed to be constant, also it can be calculated as a function of volumetric strain and other parameters can be defined by void ratio.Byrne (1991) presented an equation which correspond irrecoverable volume change and engineering shear strain with two constants.In this model, a soil mass with liquefaction potential was modelled using (N 1 ) 60 parameter as a main factor to the Finn model, so all of the soil properties needed for the model were defined for the program by (N 1 ) 60 .
Where C 1 and C 2 are constants.In many cases, C 2 = 0.4C 1 , so Eq. ( 1) involves only one independent constant.However, both C 1 and C 2 have been retained for generality (Itasca, 2008).Byrne (1991) notes that the constant, C 1 , can be derived from relative densities, as follows: = 7600 .
(2) Further, using an empirical relation between Dr and normalized standard penetration test values, (N 1 ) 60 , Then: In this study, C 2 is calculated from C 2 = 0.4C 1 .

Boundary Conditions
The seismic input is normally represented by plane waves propagating upward through the underlying material.The boundary conditions at the sides of the model must account for the free-field motion that would exist in the absence of the structure.These boundaries should be placed at distances sufficient to minimize wave reflections and achieve free-field conditions.For soils with high material damping, this condition can be obtained with a relatively small distance.However, when the material damping is low, the required distance may lead to an impractical model (Seed et al. 1975).An alternative procedure to assume that the regions far from the area of interest extended to infinity for minimizing the computation time as well as avoiding the outwards propagating waves form reflecting to the model using by free-field boundaries as shown in Fig. 2.

Element Size
Numerical distortion of propagating wave can occur in dynamic analysis as a function of modelling condition.The numerical accuracy of wave transmission is affected by both the frequency content of input wave and the wave speed characteristics of system.Kuhlemeyer & Lysmer (1973) showed that for an accurate representation of wave transmission through the soil model, the spatial element size should be smaller than 1/10 to 1/8 of the wavelength associated with the highest frequency component of input wave i.e.

∆ = 9
(5) where, λ = wave length associated with the highest frequency component that contains significant energy.Considering the above mentioned criteria, the element size is defined small enough to allow the seismic wave propagation throughout the analysis.

Damping
Material damping in soil is generally because of its viscosity, friction and plasticity development (Ebrahimian, 2012).The damping in the numerical dynamic analysis should reproduce the energy losses in the natural system when subjected to a dynamic loads.In this study, the dynamic damping is provided by Hysteretic damping combined by small amount of Rayleigh damping of 0.2%.Hysteretic damping was used as backbone curve by Seed and Idriss (1970) for granular material and Sun et al. (1988) for cohesive material.Variables included the embankment height, relative density of embankment, peak acceleration of ground motion, mechanism earthquake, and site class of base layer.The embankment was modelled as a uniform zone of cohesionless soil exhibiting a low to high relative density (specifically modeled as a sand with N = 8, 16, and 32 blows/ft) and corresponding angle of internal friction (see Table 1).The configuration of the model embankment used is shown in Fig. 5.
Each site condition was evaluated with a suite of six acceleration time histories including three different PGA and two different mechanism earthquake.
The main objective of this study was to develop a simplified approach for estimating lateral deformations and predicting liquefaction vulnerability of embankment.Accordingly, a series of parametric studies were performed numerically by using finite difference method FLAC 2D.Pertinent results consisted of excess pore pressure generation and deformations within the entire soil mass and also provided in Chapter 5.
The magnitude and direction of the deformations are computed and a contour plot of deformation can be generated for the embankment and affected base layer.The slope deformations computed using numerical modelling method must account for the fact that the numerical model provides the variation of displacement throughout the slope.In this study, the maximum computed deformation is used.

Embankment Geometry
Fig. 5 illustrates the finite difference mesh adopted in the analyses.These elements used in the mesh was about 1 x 1 m and has been accounted the maximum grid size as suggested by Kuhlmeyer (1973).Embankment heights of 5 and 10 m were modelled.The ground water level were modelled 4 meter above ground surface in the both embankment height case.Toe of the slope was modelled by rigid material consisted of rock.Width and height base layer consisted clay material were modelled by 200 and 30 meter, respectively.
In the static analyses, displacements along the vertical boundaries were restrained in horizontal direction, while displacement at the base were restricted in both the horizontal and vertical direction.In the dynamic analyses, (i) free-field boundaries were utilized on the left and right boundaries to minimize wave reflections, which otherwise would reflect back into the model, and (ii) quiet boundary were used at the bottom of the model.
In this study, the earthquake input motion must be applied as a stress boundary.Otherwise, the effect of the quiet boundary will be nullified when the input is applied as acceleration or velocity wave.

Material Properties
The material properties used in the analyses are shown in Fig. 5 and summarized in Table 1.The embankment was modelled as consisting of sand material as common material in many field applications.In this study, sand densities of 40, 60 and 80% were analysed in order to investigate the relationships among relative densities of sand, excess pore-pressures and lateral displacement of embankment.
To conduct a FLAC analysis, the soil density (γ), angle of internal friction (φ), cohesion (c), bulk modulus (K), shear modulus (G), Poisson's ratio (μ), and corrected penetration resistance ((N 1 ) 60 ) are also needed.The stiffness and strength parameters of all material were estimated from the relative densities using wellestablished correlations.

Ground Motions
In the absence of record ground motion in Indonesia, synthetic earthquake ground motion is used similar to the characteristics of Jakarta earthquake with 2475 return period.Synthetic earthquake ground motion is obtained from spectral matching to the original ground motion using the EZ-FRISK program.The original motion input is selected from several earthquake records in several locations outside Indonesia.The ground motion is referring to SSRA and Time History Generation Report of the Thamrin Nine Building, 2016.  (1Dr = Relative density of sand material (2) (N1)60 = Normalized SPT values to 1 atm and 60% blow count energy (3) γm = Bulk density (4) G & K = Shear and bulk modulus (5) C1 & C2 = Constant of finn-byrne constitutive model In this study, two different mechanism such shallow crustal and megathrust mechanisms are selected and scaled linearly to PGAs 0.1, 0.3, and 0.5 g.However, the frequency motions are scaled based on the time step (Δt) of the predominant period target.Furthermore, the scaled ground motions transferred by standard amplification analysis using NERA software and applied as stress boundary at the base of numerical model which is 30 m below the ground surface.The corresponding acceleration time histories are shown in Fig. 7 and 8 .

Numerical Results and Discussions
In this study, the effects of soil type (in terms of relative density), earthquake intensity, embankment height, and site class on liquefaction susceptibility of the soil were evaluated.The horizontal deformation and excess pore water pressure generation were also investigated.The results based on numerically simulation are discussed as follows.
The excess pore water pressure ratio was defined in the software by Fish function as shown in Equation 6, in which the soil was in the state of liquefaction for r u =1.
Where Δu is the difference between the current and hydrostatic pore pressure, and σ' vo is the initial vertical effective stress.Fig. 9 illustrate typical excess pore water pressure ratio at the end of the earthquake, while Fig. 10 describe typical time history of excess pore pressure generation.It was observed that the liquefaction susceptibility of the soil decreases in depth with increase in the soil relative density.As expected, by increasing the relative density, the excess pore pressure generation decreased.Significant decrease is shown in sand relative density of 80% to 60% and 40%.The soil deformations at sand material near the toe of the slope were found to be the highest horizontal deformation.At sand material near the toe of the slope, the horizontal displacements of the embankment increase with increasing embankment height.While at sand material far the toe of the slope, the higher of the embankment produce to smaller deformation and pore pressure.These results may be caused by the effect of the initial effective stress generation and the geometry factor of the model.
Typical horizontal deformation at the end of the earthquake and time history of horizontal deformation are shown in Fig. 11 and 12, respectively.In general, the higher sand density material would reduce the results of horizontal deformation and excess-water pressure-ratio (r u ).Sand with the value of Dr = 40% above the SE site class lead to liquefy (r u =1) in all analysis schemes at PGA of 0.1g to 0.5g.While, liquefaction does not occur at the value of Dr = 80% at PGA 0.1g to 0.5g.However, the horizontal deformation needs to be considered in the design process.
It can be observed that increase in input motion intensity leads to significant increase in displacements and its excess-pore-pressure generation.Horizontal deformations show slightly higher value in Megathrust compared to Shallow Crustal mechanism earthquake.
Coupled effects of embankment height, maximum acceleration at the base layer, soil density of embankment, and horizontal deformation & excess pore pressure ratio is presented by chart in Fig. 13 to 20.Both horizontal deformation and excess pore water pressure in relative density of 50% and 70% are estimated by average of the closest value.

Concluding Remarks
Parametric study on dynamic analysis involving 72 simulation models have been conducted adopting effective stress analyses approach.Various geometric configurations, soil material, and earthquake inputmotions have been made.
Effects on parameters such as soil density and embankment height are studied on the horizontal permanent deformations and generation of excess porewater pressure.Results of this study are presented in the form of a chart, as a practical reference in designing height of embankment and sand density for SD or SE site class commonly found in practice.This chart is useful for estimating density and embankment height in preliminary design of practical engineering process.
This chart facilitates reasonable estimate of lateral displacement and maximum pore-pressure ratio of an embankment subjected to various levels of groundmotion intensity.

Fig. 5 .
Fig. 5. Typical Geometry and Material Properties used in The Analyses

Fig. 9 .Fig. 10 .
Fig. 9. Typical Excess Pore Water Pressure Ratio Contour of 5 m Height Embankment due to 0.3g Megathrust Earthquake at Relative Density Sand of 40% above SD Site Class

Fig. 11 .Fig. 12 .
Fig. 11.Typical Horizontal Deformation Contour of 5 m Height Embankment due to 0.3g Megathrust Earthquake at Relative Density Sand of 40% above SD Site Class

Fig. 13 .
Fig. 13.Maximum Horizontal Deformation at 5 m Embankment Height above SE Site Class

Fig. 14 .
Fig. 14.Maximum Horizontal Deformation at 5 m Embankment Height above SD Site Class

Fig. 20 .
Fig. 20.Maximum Excess Pore Pressure Ratio at 10 m Embankment Height above SD Site Class

Table 1 .
Material Properties used in The Analyses