Assessment of Pohang Earthquake-Induced Liquefaction at Youngil-Man Port Using the UBCSAND2 Model

The practical constitutive model UBCSAND2, which combines two-mobilized planes—a maximum shear stress plane and a horizontal plane within a framework of classical plasticity approach—is used to incorporate shear-induced effects in both loading and unloading as well as principal stress rotation effects. UBCSAND2 was calibrated by capturing cyclic direct simple shear (CDSS) test results of Pohang sand, which was collected from liquefied paddy fields due to the 2017 Pohang earthquake (Mw = 5.4) in South Korea. The model procedure focuses on simple shear condition because it best simulates field conditions under earthquake loading. The calibrated UBCSAND2 model is then used to assess the liquefaction-induced damages that occurred at the quay wall and backfill layer in Youngil-man port near the epicenter of the Pohang earthquake. The numerical results show that liquefaction mostly occurred in silty sand layers, in which the excess pore pressure ratio reached almost one. The estimated displacements of the quay wall and the predicted settlement of reclaimed area obtained from the analysis were in good agreement with those obtained from field measurements.


Introduction
Earthquake-induced liquefaction is one of the primary causes of major damages to geotechnical structures under saturated soil conditions. Since the 1950s, very complex soil nonlinear models have been analytically developed to predict the phenomenon [1]. However, to date a reliable prediction of liquefaction-resulted displacements that developed in the plasticity formulation has still been challenging. Beaty and Perlea [2] pointed out these challenges should be considered when choosing and applying a constitutive model for a finite element analysis to predict liquefaction, which include: (i) adequate capture of cyclic soil element responses such as shear stress-strain curve, effective stress path, generation of pore pressure, (ii) best fit in the relationship between cyclic stress ratio (CSR = τ xy /σ y ) and number of cycles to liquefaction (N liq ) obtained from laboratory test results and predictions, (iii) reasonable agreement of the calibrations to the empirical relationships (i.e., field test data) for liquefaction triggering, (iv) transparent input parameters based on the realistic properties of soils (e.g., relative density D r ), and (v) robust application to many elements analysis via back analysis of case histories. Following this, the fully coupled effective stress model is considered to be the most sophisticated class of continuum models of liquefaction [2]. Another one of the conventional models that has been widely used in earthquake engineering applications (e.g., dams, bridges, buildings, and ports) is the UBCSAND model. A major advantage of this model is its relative simplicity in predicting 2.1. 1

. Elastic Response
The elastic response is assumed to be isotropic and composed from shear modulus G e and bulk modulus B e . It is also a function of mean effective stress as follows: G e = k e G ·P a · σ m P a me (2) where k e B and k e G are the elastic bulk and shear modulus number, respectively. ne and me are the stress exponent for B e and G e (approximately 0.5 for both), respectively. P a is the atmospheric pressure (approximately 100 kPa). σ m is the mean effective stress.

Yield Surface
UBCSAND2 model is formulated in the framework of classical elasto-plasticity, in which yield locus comprises the two-mobilized plane: (i) maximum shear stress plane and (ii) horizontal plane in the stress space of shear stress τ and effective normal stress σ as shown in Figure 1. The failure state line is the yield locus with the highest stress ratio that corresponds with the Mohr-Coulomb failure line (represented by a slope of sin φ f ). The line with slope of sin φ cv is the yield locus with a stress ratio that corresponds to the phase transformation from failure state to residual state of sand (represented by a slope of sin φ cv ). Yield loci for both planes f i are specified in the same way as follows: where i indicates the maximum shear plane (i = 1) or the horizontal plane (i = 2), t i is the shear stress acting on each mobilized plane (i.e., t 1 = τ = (σ x − σ y )/2, t 2 = τ xy ), σ m is the mean effective stress (i.e., (σ x + σ y )/2), and φ mi is the friction angle mobilized on each plane. In the UBCSAND2 model, both maximum shear stress and horizontal plane work at the same time for loading when only the horizontal plane is triggered upon unloading. Therefore, the modification of shear stress ratio is essential to account for the contribution of each plane to the plastic deformation as shown in Figure 2. First, the initial yield locus is located at the stress reversal point C as shown Figure 2 and the yield locus enters into unloading phase until the shear stress changes signs, or reversal occurs. At this time, the mobilized shear stress and shear stress at failure are modified to calculate the plastic shear modulus. (i.e., τ * xy = τr -τxy for unloading and τ * xy = τr + τxy for reloading and τ * f = τr + τf , where τr is the shear stress at failure). The stress ratio on each mobilized plane ηi is expressed as follows:  In the UBCSAND2 model, both maximum shear stress and horizontal plane work at the same time for loading when only the horizontal plane is triggered upon unloading. Therefore, the modification of shear stress ratio is essential to account for the contribution of each plane to the plastic deformation as shown in Figure 2. First, the initial yield locus is located at the stress reversal point C as shown Figure 2 and the yield locus enters into unloading phase until the shear stress changes signs, or reversal occurs. At this time, the mobilized shear stress and shear stress at failure are modified to calculate the plastic shear modulus. (i.e., τ * xy = τ r − τ xy for unloading and τ * xy = τ r + τ xy for reloading and τ * f = τ r + τ f , where τ r is the shear stress at failure). In the UBCSAND2 model, both maximum shear stress and horizontal plane work at the same time for loading when only the horizontal plane is triggered upon unloading. Therefore, the modification of shear stress ratio is essential to account for the contribution of each plane to the plastic deformation as shown in Figure 2. First, the initial yield locus is located at the stress reversal point C as shown Figure 2 and the yield locus enters into unloading phase until the shear stress changes signs, or reversal occurs. At this time, the mobilized shear stress and shear stress at failure are modified to calculate the plastic shear modulus. (i.e., τ * xy = τr -τxy for unloading and τ * xy = τr + τxy for reloading and τ * f = τr + τf , where τr is the shear stress at failure). The stress ratio on each mobilized plane ηi is expressed as follows:  The stress ratio on each mobilized plane η i is expressed as follows: When the yield locus moves from A to B in Figure 1, it induces a plastic shear strain increment dγ p 1 related to the shear stress ratio η 1 . It is assumed that the relationship of stress ratio η 1 and plastic strain γ p 1 is a hyperbolic curve, which is controlled by plastic shear modulus G p 1 as shown in Figure 3. This correlation can be expressed as follows: where dγ p 1 is the plastic shear strain increment, G p 1 is the plastic shear modulus, dη 1 is the stress ratio increment.
where dγ p 1 is the plastic shear strain increment, G p 1 is the plastic shear modulus, dη1 is the stress ratio increment.
G p 1 is the plastic shear modulus and given by where η1 is the stress ratio at failure (= sin ϕf), Rf is the failure ratio for the hyperbolic relationship (= ηf/ηult). Rf varies from 0.6 to 1.0 and decreases with increasing relative density. G p max is the plastic modulus at η = 0 and can be expressed as follows: where k p G is the plastic shear modulus number. np is the stress exponent for G p max (approximately 0.4).  Herein, the position of the yield locus dϕ m1 can be calculated as follows: The flow rule acting on maximum shear stress plane, which is obtained from the well-known Rowe stress-dilatancy theory [21], is given as follows: G p 1 is the plastic shear modulus and given by where η 1 is the stress ratio at failure (= sin φ f ), R f is the failure ratio for the hyperbolic relationship (= η f /η ult ). R f varies from 0.6 to 1.0 and decreases with increasing relative density. G p max is the plastic modulus at η = 0 and can be expressed as follows: where k p G is the plastic shear modulus number. np is the stress exponent for G p max (approximately 0.4). Herein, the position of the yield locus dφ m1 can be calculated as follows: The flow rule acting on maximum shear stress plane, which is obtained from the well-known Rowe stress-dilatancy theory [21], is given as follows: where dε v1 p is the plastic volumetric strain increment, φ cv is the friction angle at constant volume, η 1 is the current stress ratio, and ψ 1 is the dilatancy angle.

Hardening Rule and Flow Rule on the Horizontal Plane
For loading phase, the hardening rule and flow rule are similar for both maximum shear stress and horizontal planes. The only difference is the stress ratio η.
Upon loading, plastic deformation is controlled by conditions on the horizontal plane, in which consistence condition of its yield surface is derived from Equation (3) as follows: The plastic shear modulus G p 2 of the horizontal plane for unloading phase is expressed as follows: where η * and η f * are modified stress ratios based on modified shear stresses τ * xy and τ * f as shown in Figure 2.

Hardening of Stress-Dilatancy
In order to account for reducing volumetric strains with number of cycles (volumetric hardening), it is necessary to either change the stress-dilatancy equation (e.g., state-dependent dilatancy) or stiffen the plastic shear modulus. Hence, the stress-dilatancy relationship in Equation (9) is modified by the multiplication to a cyclic hardener as the following expression: where C h is a cyclic hardening parameter that varies with relative density. From experimental data of Silver and Seed [22] on crystal silica sand, and Sriskandakumar [23] on Fraser river sand, the value of C h can be expressed as in Equation (14). Additionally, γ is the current shear strain.

UBCSAND2 Calibration for Undrained CDSS Tests on Pohang Sand
After the Pohang earthquake, many sand boils occurred in the paddy fields near the epicenter. The ejected sands from these sand boils were collected and tested by a series of undrained CDSS tests [24]. These results were used to calibrate the element response of the UBCSAND2 model. The samples were reconstituted with an air-pluviation method into loose (D r = 40%) and dense (D r = 80%) states under vertical effective stresses of 100 and 200 kPa, respectively. Symmetrical and non-symmetrical loading (i.e., α = τ sta /σ' y0 = 0, and α > 0) were applied in the CDSS tests to investigate the effects of initial shear stress on the cyclic resistance of sand. Figure 4 shows the experimental results for the relationship of cyclic stress ratio (CSR) and number of cycles to liquefaction (N liq ) of Pohang sand in both loose and dense states under a vertical effective stress (σ v0 ) of 100 kPa. At N liq = 15, if the relative density of sand increased from 40% to 80%, then the resistance of liquefaction increased by approximately 28.6% (from 0.14 to 0.18). Based on the SPT base clean-sand curve from a workshop of the National Center for Earthquake Research (NCEER) summarized by Youd et al. [25], the corrected SPT values N 1(60) of 12.5 and 16.5 were inferred from a CSR of 0.14 and 0.18 at σ v0 = 100 kPa, respectively.  From Figure 4, two corrected SPT blow counts N1(60) = 12.5 for loose sand and N1(60) = 16.5 for dense sand were used to calibrate the UBCSAND2 model. The generic input parameters of the UBCSAND2 model are the same as those that were used in the commercial version v.904aR [3], which is shown in Table 1. We note that the UBCSAND2 model has an additional input parameter corresponding to the cyclic hardening parameter Ch, which is a function of relative density Dr based on the experiment results of Fraser river sand. A detailed description of this parameter was previously presented. The calibration parameters were calculated identically to the generic parameters for both loose and dense sand, to match the laboratory data. It is interesting that providing only corrected SPT blow counts N1(60) as the input parameters of UBCSAND2, and the bulk modulus of water, a predictive model of liquefaction can be achieved. Figures 5 and 6 show the comparison results between the laboratory tests and the calibrated UBCSAND2 model for loose and dense Pohang sand with σ´v0 = 100 kPa, CSR = 0.14, and α = 0 (i.e., without initial static shear stress), respectively. The predictions are generally in good agreement with the measurements in terms of effective stress path (i.e., vertical effective stress versus shear stress), shear stress-strain curve, pore water pressure generation, and shear strain versus number of cycles. In the effective stress path, there are two distinct phases: (i) rapid effective drop for the initial three cycles of loose sand and five cycles of dense sand, and (ii) slow effective stress drop for 10 cycles of loose sand and 32 cycles of dense sand. This leads to a rapid rise in pore pressure of up to 35-45% the initial vertical effective stress during three cycles and five cycles for loose and dense sand, respectively. Herein, the effect of PSR is accounted for in the model, in which both the maximum shear stress plane and horizontal plane participate in the generation of plastic shear strains. This effect gradually reduces as the direction of the principal stress rotates, and it vanishes when σ'x = σ'y (i.e., at three cycles for loose sand and four cycles for dense sand). In other words, during phase (i), the two planes work together due to PSR, and then only the maximum plane works in phase (ii) for the remaining cycles because no more PSR occurs.  From Figure 4, two corrected SPT blow counts N 1(60) = 12.5 for loose sand and N 1(60) = 16.5 for dense sand were used to calibrate the UBCSAND2 model. The generic input parameters of the UBCSAND2 model are the same as those that were used in the commercial version v.904aR [3], which is shown in Table 1. We note that the UBCSAND2 model has an additional input parameter corresponding to the cyclic hardening parameter C h , which is a function of relative density D r based on the experiment results of Fraser river sand. A detailed description of this parameter was previously presented. The calibration parameters were calculated identically to the generic parameters for both loose and dense sand, to match the laboratory data. It is interesting that providing only corrected SPT blow counts N 1(60) as the input parameters of UBCSAND2, and the bulk modulus of water, a predictive model of liquefaction can be achieved. Figures 5 and 6 show the comparison results between the laboratory tests and the calibrated UBCSAND2 model for loose and dense Pohang sand with σ v0 = 100 kPa, CSR = 0.14, and α = 0 (i.e., without initial static shear stress), respectively. The predictions are generally in good agreement with the measurements in terms of effective stress path (i.e., vertical effective stress versus shear stress), shear stress-strain curve, pore water pressure generation, and shear strain versus number of cycles. In the effective stress path, there are two distinct phases: (i) rapid effective drop for the initial three cycles of loose sand and five cycles of dense sand, and (ii) slow effective stress drop for 10 cycles of loose sand and 32 cycles of dense sand. This leads to a rapid rise in pore pressure of up to 35-45% the initial vertical effective stress during three cycles and five cycles for loose and dense sand, respectively. Herein, the effect of PSR is accounted for in the model, in which both the maximum shear stress plane and horizontal plane participate in the generation of plastic shear strains. This effect gradually reduces as the direction of the principal stress rotates, and it vanishes when σ x = σ y (i.e., at three cycles for loose sand and four cycles for dense sand). In other words, during phase (i), the two planes work together due to PSR, and then only the maximum plane works in phase (ii) for the remaining cycles because no more PSR occurs.

Figure 5.
Comparison of the UBCSAND2 model for loose Pohang sand and laboratory CDSS test at CSR = 0.14, σ v0 = 100 kPa, and α = 0.0 in the variations of (a) stress path, (b) shear stress-strain curve, (c) excess pore water pressure generation, (d) shear strain and number of cycles to liquefaction.
Appl. Sci. 2020, 10, x FOR PEER REVIEW 9 of 26 Figure 5. Comparison of the UBCSAND2 model for loose Pohang sand and laboratory CDSS test at CSR = 0.14, σ´v0 = 100 kPa, and α = 0.0 in the variations of (a) stress path, (b) shear stress-strain curve, (c) excess pore water pressure generation, (d) shear strain and number of cycles to liquefaction. The typical comparisons between the undrained CDSS test results and calibrated UBCSAND2 model considering the initial static shear stress bias (i.e., α = τsta/σ'y0 > 0) are shown in Figures 7 and 8. The level of the initial static shear stress ratio α is defined as the ratio of the driving initial static shear stress, τsta, to the initial vertical effective stress, σ'y0 (i.e., α = τsta/σ'y0). This kind of effect is varied in the soil elements beneath the slope. The predictions generally give a reasonable representation for stress reversal patterns, including an accumulating pattern of shear strain. The typical comparisons between the undrained CDSS test results and calibrated UBCSAND2 model considering the initial static shear stress bias (i.e., α = τ sta /σ y0 > 0) are shown in Figures 7 and 8. The level of the initial static shear stress ratio α is defined as the ratio of the driving initial static shear stress, τ sta , to the initial vertical effective stress, σ y0 (i.e., α = τ sta /σ y0 ). This kind of effect is varied in the soil elements beneath the slope. The predictions generally give a reasonable representation for stress reversal patterns, including an accumulating pattern of shear strain. Figure 9 compares the CSR-N liq curves from the CDSS tests and UBCSAND2-calibrated results for Pohang sand in loose and dense states with σ v0 = 100 and 200 kPa, respectively. There is a good agreement between the trend line from the UBCSAND2-clibrated results and CDSS test results. The effect of confining stress on CSR is also shown, in which the increase of the vertical effective stress from 100 to 200 kPa results in the decrease of CSR for both density states of Pohang sand. However, this effect is more significant in loose sand compared to that in dense sand.
From the calibration results for the undrained CDSS tests on Pohang sand, the generic input parameters of the UBCSAND2 model depicted in Table 1 are used to numerically assess the liquefaction from the Pohang earthquake in the next section.     effect of confining stress on CSR is also shown, in which the increase of the vertical effective stress from 100 to 200 kPa results in the decrease of CSR for both density states of Pohang sand. However, this effect is more significant in loose sand compared to that in dense sand. From the calibration results for the undrained CDSS tests on Pohang sand, the generic input parameters of the UBCSAND2 model depicted in Table 1 are used to numerically assess the liquefaction from the Pohang earthquake in the next section.

Study Area and Ground Motion of the Pohang Earthquake
The Pohang earthquake occurred on 15 November 2017, causing damage to the caisson-type quay wall structures and reclaimed area at the Youngil-man port, which is located about 6 km from the epicenter, as shown in Figure 10. After the earthquake, Kim et al. [11] conducted a field study and indicated that the settlement of reclaimed area was approximately 10-20 cm, the caisson horizontal displacement was in a range of 5-15 cm, and the caisson settlement was less than 10 cm. Additionally, many sand boils could be found in this area, as shown in Figure 10.
During the Pohang earthquake, ground motions were recorded at four stations-PHA2, CHS, DKJ, and HAK-with distances from the epicenter of 9, 26, 29, and 23 km, respectively, as shown in Figure 10. The seismic motions were obtained from the National Earthquake Comprehensive Information Center [26]. Table 2 summarizes the information from the four selected stations: peak ground acceleration (PGA)-north-south (NS), east-west (EW), and up-down (UD)-locations, distances from epicenter, and soil classifications. The site condition in Pohang is classified as SB based on the Korea Building Code [27]. This means that an average shear wave velocity for the ground of up to a depth of 30 m (i.e., Vs 30 ) is 760-1500 m/s. We note that the horizontal peak accelerations in the EW and UD directions were significantly lower than those recorded in the NS direction for all of stations. Furthermore, the peak accelerations decreased remarkably when the stations' distances from the epicenter increased.

Study Area and Ground Motion of the Pohang Earthquake
The Pohang earthquake occurred on 15 November 2017, causing damage to the caisson-type quay wall structures and reclaimed area at the Youngil-man port, which is located about 6 km from the epicenter, as shown in Figure 10. After the earthquake, Kim et al. [11] conducted a field study and indicated that the settlement of reclaimed area was approximately 10-20 cm, the caisson horizontal displacement was in a range of 5-15 cm, and the caisson settlement was less than 10 cm. Additionally, many sand boils could be found in this area, as shown in Figure 10.   During the Pohang earthquake, ground motions were recorded at four stations-PHA2, CHS, DKJ, and HAK-with distances from the epicenter of 9, 26, 29, and 23 km, respectively, as shown in Figure 10. The seismic motions were obtained from the National Earthquake Comprehensive Information Center [26]. Table 2 summarizes the information from the four selected stations: peak ground acceleration (PGA)-north-south (NS), east-west (EW), and up-down (UD)-locations, distances from epicenter, and soil classifications. The site condition in Pohang is classified as S B based on the Korea Building Code [27]. This means that an average shear wave velocity for the ground of up to a depth of 30 m (i.e., V s 30 ) is 760-1500 m/s. We note that the horizontal peak accelerations in the EW and UD directions were significantly lower than those recorded in the NS direction for all of stations. Furthermore, the peak accelerations decreased remarkably when the stations' distances from the epicenter increased. Note: E = east, N = north, NS = north-south, EW = east-west, UD = up-down, PGA = peak ground acceleration.
From Figure 10, the observed site (i.e., Youngil-man port) locates nearly 6 km from the epicenter, which is closer than all the recording stations. However, there is no seismic instruments available in this area because they are only installed in the rock formation that is visible on the surface (i.e., outcrop). If near-field effects are considered, the ground motions should be modified based on distances from the epicenter to the selected station, and the observed site. This kind of modification requires a modern and complex model based on fault rupture surfaces, geometries, and energy radiation distributions. For example, the Southern California Earthquake Center Broadband Platform (SCEC-BPP) was successfully used to address the near-field effects of the 1971 San Fernando earthquake-induced damages in Lower and Upper San Fernando dams [28,29]. For the Pohang earthquake, neither reliable research nor available data have been performed to study the near-field effects to the recorded PGAs. Therefore, this study applied a simple method to scale the PGAs from the available data provided by the National Earthquake Comprehensive Information Center [26]. A relationship between recorded PGAs and distances from the epicenter of different recording stations was established in Figure 11. Based on the logarithm curve, the PGA of 0.350 g at Youngil-man port was obtained. Then, the acceleration-time history recorded at PHA2 station was scaled up to 0.35 g and it is considered as the target earthquake given as outcrop motion at the observed site.
For FLAC analyses [30], the seismic input must be applied at the base of the model while the target earthquake at Youngil-man port is the outcrop motion. The appropriate input motion at depth can be computed through a deconvolution analysis using the one-dimensional (1-D) wave propagation code such as the equivalent-linear program PROSHAKE [31] as shown in Figure 12. A typical soil profile at the observed site up to 50 m of depth is presented in Figure 12a, in which the data is obtained from Kim et al. [11] and Lee et al. [12]. Under a depth of 30 m, the bedrock is primarily distributed up to 500 m of depth as the report of Lee et al. [32] and Kang et al. [9]. Additionally, the corrected N 1(60),cs is calculated from the data provided by the Integrated DB Center of National Geotechnical Information [33]. To simply the input parameters for both PROSHAKE and FLAC analysis (i.e., N 1(60),cs ), Kim et al. [11] and Lee et al. [12] calculated the corresponding N 1(60),cs values of sandy soil layers by the experimental works of Mun [34] and Park et al. [35] in relation to the liquefied soils due to the Pohang earthquake. The calculation procedure for the N 1(60),cs value was based on a method proposed by Mortia et al. [36] and Iai et al. [37]. The sediment and weather rock layer are assumed as N 1(60),cs = 20 and N 1(60),cs = 30, respectively. Although the bedrock is very hard for the SPT test, it can be assumed as N 1(60),cs = 50 for the equivalent-linear analysis. From Figure 12b, there is a reasonable agreement between the corrected N 1(60),cs from SPT test and the simplification. Based on the simplified N 1(60),cs , the maximum shear modulus G max can be determined by the empirical function proposed by Seed et al. [38] as in Equations (15) and (16). If shear modulus G max and unit weight γ are inputted, shear wave velocity V s would be automatically calculated by PROSHAKE software. It can be seen that V s is varied in the ranges 180-270 m/s for silty sand, 285 m/s for sediment, 300-330 m/s for weather rock, and 760 m/s for bedrock layer.
where G max is maximum shear modulus, K 2max is shear modulus number, P a is atmospheric pressure (≈100 kPa), σ m is mean effective stress, N 1(60),cs is corrected SPT number of blow counts.
Appl. Sci. 2020, 10, x FOR PEER REVIEW 13 of 26 Figure 11. Estimated peak ground acceleration (PGA) at Youngil-man port from the relationship between recorded PGAs and distances from the epicenter of different recording stations.
For FLAC analyses [30], the seismic input must be applied at the base of the model while the target earthquake at Youngil-man port is the outcrop motion. The appropriate input motion at depth can be computed through a deconvolution analysis using the one-dimensional (1-D) wave propagation code such as the equivalent-linear program PROSHAKE [31] as shown in Figure 12. A typical soil profile at the observed site up to 50 m of depth is presented in Figure 12a, in which the data is obtained from Kim et al. [11] and Lee et al. [12]. Under a depth of 30 m, the bedrock is primarily distributed up to 500 m of depth as the report of Lee et al. [32] and Kang et al. [9]. Additionally, the corrected N1(60),cs is calculated from the data provided by the Integrated DB Center of National Geotechnical Information [33]. To simply the input parameters for both PROSHAKE and FLAC analysis (i.e., N1(60),cs), Kim et al. [11] and Lee et al. [12] calculated the corresponding N1(60),cs values of sandy soil layers by the experimental works of Mun [34] and Park et al. [35] in relation to the liquefied soils due to the Pohang earthquake. The calculation procedure for the N1(60),cs value was based on a method proposed by Mortia et al. [36] and Iai et al. [37]. The sediment and weather rock layer are assumed as N1(60),cs = 20 and N1(60),cs = 30, respectively. Although the bedrock is very hard for the SPT test, it can be assumed as N1(60),cs = 50 for the equivalent-linear analysis. From Figure 12b, there is a reasonable agreement between the corrected N1(60),cs from SPT test and the simplification. Based on the simplified N1(60),cs, the maximum shear modulus Gmax can be determined by the empirical function proposed by Seed et al. [38] as in Equations (15) and (16). If shear modulus Gmax and unit weight γ are inputted, shear wave velocity Vs would be automatically calculated by PROSHAKE software. It can be seen that Vs is varied in the ranges 180-270 m/s for silty sand, 285 m/s for sediment, 300-330 m/s for weather rock, and 760 m/s for bedrock layer.
where Gmax is maximum shear modulus, K2max is shear modulus number, Pa is atmospheric pressure Figure 11. Estimated peak ground acceleration (PGA) at Youngil-man port from the relationship between recorded PGAs and distances from the epicenter of different recording stations. Figure 13 presents the numerical results obtained from the 1-D equivalent-linear analysis in PROSHAKE. By the application of deconvolution procedure as shown in Figure 13a, the outcrop motion at surface is inputted at the base at a depth of 50 m. Then, the wave propagation would be generated to the above layers depending on their properties. Three ground motions at the depth of 50 m, 25 m, and free-surface are extracted for the comparison as shown in Figure 13b-d. Although the PGA at the free-surface is nearly the same as PGA at the base, however, there is a significantly higher acceleration pulse at the ground surface. At the depth of 25 m, the maximum acceleration is approximately 0.252 g, in which the seismic energy applied at the base reduced almost 30%.
The extracted ground motion at the depth of 25 m is selected as the input ground motion for FLAC analysis in the next section. Figure 14 depicts the acceleration-time history, velocity-time history, displacement-time history, and power spectrum of the acceleration of the Pohang earthquake, respectively. The acceleration was extracted to approximately 40 s, filtered to 20 Hz, and corrected by a baseline drift. Then, it was implemented in the numerical analysis to assess the Pohang earthquake-induced liquefaction at Youngil-man port.   Figure 13a, the outcrop motion at surface is inputted at the base at a depth of 50 m. Then, the wave propagation would be generated to the above layers depending on their properties. Three ground motions at the depth of 50 m, 25 m, and free-surface are extracted for the comparison as shown in Figure 13b-d. Although the PGA at the free-surface is nearly the same as PGA at the base, however, there is a significantly higher acceleration pulse at the ground surface. At the depth of 25 m, the maximum acceleration is approximately 0.252 g, in which the seismic energy applied at the base reduced almost 30%. The extracted ground motion at the depth of 25 m is selected as the input ground motion for FLAC analysis in the next section. Figure 14 depicts the acceleration-time history, velocity-time history, displacement-time history, and power spectrum of the acceleration of the Pohang earthquake, respectively. The acceleration was extracted to approximately 40 s, filtered to 20 Hz, and corrected by a baseline drift. Then, it was implemented in the numerical analysis to assess the Pohang earthquake-induced liquefaction at Youngil-man port. The extracted ground motion at the depth of 25 m is selected as the input ground motion for FLAC analysis in the next section. Figure 14 depicts the acceleration-time history, velocity-time history, displacement-time history, and power spectrum of the acceleration of the Pohang earthquake, respectively. The acceleration was extracted to approximately 40 s, filtered to 20 Hz, and corrected by a baseline drift. Then, it was implemented in the numerical analysis to assess the Pohang earthquake-induced liquefaction at Youngil-man port.  Figure 15 presents the representative cross section of Youngil-man port, which is in the northsouth direction and was created using the FLAC2D program [30]. It was mostly built from gravel as a backfill material near the quay wall, and silty sandy soil for the reclaimed area. These layers are extremely sensitive to liquefaction under earthquake motion. Lee et al. [32] and Kang et al. [9] Figure 15 presents the representative cross section of Youngil-man port, which is in the north-south direction and was created using the FLAC2D program [30]. It was mostly built from gravel as a backfill material near the quay wall, and silty sandy soil for the reclaimed area. These layers are extremely sensitive to liquefaction under earthquake motion. Lee et al. [32] and Kang et al. [9] indicated that the bedrock of the Pohang basin, including the observed area, is approximately 500 m deep, and the mudstone is primarily distributed below a depth of 20 m up to the bedrock. In the case of Youngil-man port, the foundation soils include a thin thickness of sedimentary layer mainly from silty-like soil, and weathered rock is below the sedimentary layer [12]. In addition, a row of steel pipe piles, with a diameter of 600 mm, was installed at the depth of the sedimentary layers in order to reduce displacements if a failure occurs at the quay wall. Nodes N1, N2, N3, and N4 in Figure 15 depict the positions of the output nodes for the displacement predictions of the caisson, landfill ground behind the quay wall, and wave transmission checking. Similarly, elements E1, E2, E3, and E4 are used for the extraction of pore pressure generation and reduction in effective stress when liquefaction occurs.

Numerical Modeling
Appl. Sci. 2020, 10, x FOR PEER REVIEW 16 of 26 indicated that the bedrock of the Pohang basin, including the observed area, is approximately 500 m deep, and the mudstone is primarily distributed below a depth of 20 m up to the bedrock. In the case of Youngil-man port, the foundation soils include a thin thickness of sedimentary layer mainly from silty-like soil, and weathered rock is below the sedimentary layer [12]. In addition, a row of steel pipe piles, with a diameter of 600 mm, was installed at the depth of the sedimentary layers in order to reduce displacements if a failure occurs at the quay wall. Nodes N1, N2, N3, and N4 in Figure 15 depict the positions of the output nodes for the displacement predictions of the caisson, landfill ground behind the quay wall, and wave transmission checking. Similarly, elements E1, E2, E3, and E4 are used for the extraction of pore pressure generation and reduction in effective stress when liquefaction occurs. Material properties: Table 3 shows the input parameters implemented in the numerical analysis. In this study, two silty sand layers in the reclaimed region were considered as the liquefiable soil layers, in which a reasonable N1(60) value of 15 was determined to be the best input into the numerical analysis using the UBCSAND2 model. All of the generic parameters required in the UBCSAND2 model in Table 1 can be calculated directly via the single N1(60) value. The perfectly plastic Mohr-Coulomb model was assigned for non-liquefiable soil layers such as sediment, gravel, weather rock, and riprap. The elastic materials of concrete and steel were selected for the caisson-type quay wall Material properties: Table 3 shows the input parameters implemented in the numerical analysis. In this study, two silty sand layers in the reclaimed region were considered as the liquefiable soil layers, in which a reasonable N 1(60) value of 15 was determined to be the best input into the numerical analysis using the UBCSAND2 model. All of the generic parameters required in the UBCSAND2 model in Table 1 can be calculated directly via the single N 1(60) value. The perfectly plastic Mohr-Coulomb model was assigned for non-liquefiable soil layers such as sediment, gravel, weather rock, and riprap. The elastic materials of concrete and steel were selected for the caisson-type quay wall and the steel pipe pile, respectively.  Element size: The mesh size of the analysis was selected to be ∆l = 1 m in order to ensure accurate wave transmission, and is based on the material properties in Table 3. Silty sand 1,2 has the lowest shear-wave speed C s , for an elastic shear modulus k G e and a density ρ, which can be calculated as follows: Then, the maximum frequency transmitted in the largest zone in FLAC must be higher than the maximum frequency of the input motion [30].
From Equation (18), it can be seen that the selected mesh size of 1 × 1 m can transfer the seismic wave in an accurate way when the earthquake motion is applied at the base of the numerical analysis.
Boundary conditions for static analysis: During the static analysis, both the horizontal and vertical directions were fixed against the movement at the base of the model. Two sides of the model allowed for vertical movements, but were laterally fixed to prevent horizontal movements. At the bottom, a free pore pressure condition was assigned to allow for changes in pore pressure, but it restricted flow across the boundary. The left and right sides applied pore pressure below the level of −3.00 m relative to the ground surface, which was the maximum sea level rise. Steady-state seepage analysis was then performed to develop the initial phreatic surface.
Boundary conditions for seismic analysis: A quiet boundary condition in the xand y-directions was applied to the bottom of the model to minimize the reflection of outward propagating waves back into the model and does not allow the necessary energy radiation [30,39,40]. Along the side boundaries of the grid should not be assigned as the quiet boundaries because the wave energy can leak out of the sides. Alternatively, free-field boundaries can be used in this situation [30]. During the seismic analysis, the side boundary zones were assigned to be elastic elements with appropriate stiffness (three columns of elements on each side) to prevent free-field boundary conditions from element distortions and increase the spurious reflections off the boundary. This kind of technique can minimize wave reflections if dynamic loading is applied to the model boundaries [30]. To achieve dynamic loading at the base, the ground motions were transformed to equivalent shear waves (i.e., horizontal motion) for the weathered rock layer as follows: where σ s is the applied shear stress, ρ is the mass density, C s is the velocity of the s-wave propagation through medium, and ν s is the input velocity. Viscous damping: Rayleigh damping was used for the seismic analysis in FLAC. The predominant frequency of the ground motions (= 1.05 Hz, as in Figure 14d) was considered as the center of frequency in the analysis. The numerical results obtained from the equivalent linear analysis using PROSHAKE software were extracted for silty sand, sediment, weather rock, and gravel as shown in Figure 16. The selected layer corresponds to the maximum shear strain, which is approximately 0.03% for silty sand, 0.01% for sediment, 0.01% for weather rock, and 0.01% for gravel. Then, the damping ratio and modulus reduction ratio (G/G max ) are inferred from the curves in Figure 16 for each soil type with the equivalent uniform strain (about 50% of maximum shear strain) [30]. As a result, the damping ratio is taken as 5% for silty sand, 4% for sediment, 2% for weather rock, and 5% for gravel. Whereas, the modulus reduction factor is selected as 0.7 for silty sand, 0.78 for sediment, 0.9 for weather rock, and 0.58 for gravel.
Appl. Sci. 2020, 10, x FOR PEER REVIEW 18 of 26 16 for each soil type with the equivalent uniform strain (about 50% of maximum shear strain) [30]. As a result, the damping ratio is taken as 5% for silty sand, 4% for sediment, 2% for weather rock, and 5% for gravel. Whereas, the modulus reduction factor is selected as 0.7 for silty sand, 0.78 for sediment, 0.9 for weather rock, and 0.58 for gravel.  Figure 17 shows the stress state just before the start of the earthquake. In the reclaimed region behind the steel pipe piles, the initial horizontal stress and initial shear stress ratios are mostly around 0.5 (Jaky's lateral earth pressure coefficient K0 = 1 − sin ϕ´ ≈ 0.5 for sand [41]) and 0.0, respectively. Near the quay wall, K0 is in the range of 0.2-0.4, and α ranges from 0.05 to approximately 0.1. Although there are higher values of K0 ≈ 1.0 and α ≈ 0.15 at other locations, such as at the base of the quay wall (i.e., riprap), soils in these areas are considered to be non-liquefiable materials. Thus, their contributions to liquefaction are minimal. The predictions for the static state are reasonable for the realistic conditions of the Youngil-man port before the earthquake occurred.  Figure 17 shows the stress state just before the start of the earthquake. In the reclaimed region behind the steel pipe piles, the initial horizontal stress and initial shear stress ratios are mostly around 0.5 (Jaky's lateral earth pressure coefficient K 0 = 1 − sin φ ≈ 0.5 for sand [41]) and 0.0, respectively. Near the quay wall, K 0 is in the range of 0.2-0.4, and α ranges from 0.05 to approximately 0.1. Although there are higher values of K 0 ≈ 1.0 and α ≈ 0.15 at other locations, such as at the base of the quay wall (i.e., riprap), soils in these areas are considered to be non-liquefiable materials. Thus, their contributions to liquefaction are minimal. The predictions for the static state are reasonable for the realistic conditions of the Youngil-man port before the earthquake occurred.

Pre-Earthquake Stress Conditions
behind the steel pipe piles, the initial horizontal stress and initial shear stress ratios are mostly around 0.5 (Jaky's lateral earth pressure coefficient K0 = 1 − sin ϕ´ ≈ 0.5 for sand [41]) and 0.0, respectively. Near the quay wall, K0 is in the range of 0.2-0.4, and α ranges from 0.05 to approximately 0.1. Although there are higher values of K0 ≈ 1.0 and α ≈ 0.15 at other locations, such as at the base of the quay wall (i.e., riprap), soils in these areas are considered to be non-liquefiable materials. Thus, their contributions to liquefaction are minimal. The predictions for the static state are reasonable for the realistic conditions of the Youngil-man port before the earthquake occurred.  Figure 18 presents the maximum of the excess pore water pressure ratio Ru,seis contours at the end of the earthquake. This value is the ratio of the pore water pressure seismically generated and the initial vertical effective stress, which differs from the commonly used ratio of the current pore pressure and total vertical effective stress (ru = u/σy). Beaty and Byrne [3] suggested that the value of Ru,seis can be determined as follows:

Excess Pore Water Pressure
where σ´y is the current vertical effective stress, σ´y,i is the initial (pre-earthquake) vertical effective stress, and Δucyclic is the cyclically induced change in pore pressure.
Based on the contours of Ru,seis in Figure 18, the soils are mainly liquefied in the second silty sand layer, which is more significant in the far field and near the quay wall. Figure 18. Maximum excess pore water pressure ratio (Ru,seis = 1.0 − σ´y/σ´y,i) contours at the end of shaking (dimension unit in meters). Figure 19 presents the excess pore water pressure generation during the earthquake of observed elements in Figure 15. Elements E1, E2, and E3 are defined as the liquefied second silty sand layer, Figure 17. Initial stress state at the start of the earthquake. (a) Initial horizontal stress ratio (K 0 = σ x /σ y ) contours (dimension unit in meters); (b) initial static shear stress ratio (α = τ xy /σ y ) contours (dimension unit in meters). Figure 18 presents the maximum of the excess pore water pressure ratio R u,seis contours at the end of the earthquake. This value is the ratio of the pore water pressure seismically generated and the initial vertical effective stress, which differs from the commonly used ratio of the current pore pressure and total vertical effective stress (r u = u/σ y ). Beaty and Byrne [3] suggested that the value of R u , seis can be determined as follows:

Excess Pore Water Pressure
where σ y is the current vertical effective stress, σ y,i is the initial (pre-earthquake) vertical effective stress, and ∆u cyclic is the cyclically induced change in pore pressure. Based on the contours of R u,seis in Figure 18, the soils are mainly liquefied in the second silty sand layer, which is more significant in the far field and near the quay wall. Figure 19 presents the excess pore water pressure generation during the earthquake of observed elements in Figure 15. Elements E1, E2, and E3 are defined as the liquefied second silty sand layer, and they are assigned by the UBCSAND2 model. Element E4 is considered as non-liquefied soil and is assigned by the Mohr-Coulomb model. It is clear that elements E1 and E2 reach a peak value of R u,seis ≈ 1.0, in which liquefaction occurs in these elements. Although element E3 firstly generated to the peak value of R u,seis = 0.9, however, it decreased slightly to nearly 0.6 at the end of the shaking. In contrast, the peak R u,seis of element E4 is approximately 0.4, in which no liquefaction occurs in this element.
where σ´y is the current vertical effective stress, σ´y,i is the initial (pre-earthquake) vertical effective stress, and Δucyclic is the cyclically induced change in pore pressure.
Based on the contours of Ru,seis in Figure 18, the soils are mainly liquefied in the second silty sand layer, which is more significant in the far field and near the quay wall.  Figure 19 presents the excess pore water pressure generation during the earthquake of observed elements in Figure 15. Elements E1, E2, and E3 are defined as the liquefied second silty sand layer, and they are assigned by the UBCSAND2 model. Element E4 is considered as non-liquefied soil and is assigned by the Mohr-Coulomb model. It is clear that elements E1 and E2 reach a peak value of Ru,seis ≈ 1.0, in which liquefaction occurs in these elements. Although element E3 firstly generated to the peak value of Ru,seis = 0.9, however, it decreased slightly to nearly 0.6 at the end of the shaking. In . Figure 19. Excess pore pressure ratio generated during shaking (Ru,seis). Figure 20 shows the deformed mesh of the analysis at the end of shaking. An enlargement around the quay wall is scaled up by five times to easily identify displacements. It can be seen that the top of the quay wall moved horizontally approximately 20 cm toward the ocean and settled nearly 3 cm from its starting point, which is in agreement with the observations where the wall displacements are reported to be 5-15 cm horizontally and less than 10 cm vertically [11]. It follows from this figure that the deformation of the foundation (i.e., riprap layer) leads to displacements of the quay wall. The slippage in the interface beneath the quay wall may not be the reason for the quay wall displacements, and it is mainly due to deformation in the foundation. A detail displacement description of the observed nodes (set in Figure 15) is shown in Figure 21. It can be seen that the horizontal displacements of nodes N1 and N2 are quite similar, which are around 20 cm. However, the settlement of node N2 is slightly higher than that of node N1. The settlement predictions of nodes N3 and N4 are approximately 8 and 4 cm, respectively. When nodes N1, N2, and N3 move horizontally toward the ocean (approximately 20 cm for nodes N1 and N2 and 8 cm for node N3), node N4 deforms to the right side of the model by around 2 cm. It is noticed that the displacements of soil behind the steel pipe piles (represented by node N4) are reduced compared to those of the soil in front of the row of steel piles (represented by nodes N1, N2, and N3).  Figure 20 shows the deformed mesh of the analysis at the end of shaking. An enlargement around the quay wall is scaled up by five times to easily identify displacements. It can be seen that the top of the quay wall moved horizontally approximately 20 cm toward the ocean and settled nearly 3 cm from its starting point, which is in agreement with the observations where the wall displacements are reported to be 5-15 cm horizontally and less than 10 cm vertically [11]. It follows from this figure that the deformation of the foundation (i.e., riprap layer) leads to displacements of the quay wall. The slippage in the interface beneath the quay wall may not be the reason for the quay wall displacements, and it is mainly due to deformation in the foundation. A detail displacement description of the observed nodes (set in Figure 15) is shown in Figure 21. It can be seen that the horizontal displacements of nodes N1 and N2 are quite similar, which are around 20 cm. However, the settlement of node N2 is slightly higher than that of node N1. The settlement predictions of nodes N3 and N4 are approximately 8 and 4 cm, respectively. When nodes N1, N2, and N3 move horizontally toward the ocean (approximately 20 cm for nodes N1 and N2 and 8 cm for node N3), node N4 deforms to the right side of the model by around 2 cm. It is noticed that the displacements of soil behind the steel pipe piles (represented by node N4) are reduced compared to those of the soil in front of the row of steel piles (represented by nodes N1, N2, and N3).

Displacements
horizontal displacements of nodes N1 and N2 are quite similar, which are around 20 cm. However, the settlement of node N2 is slightly higher than that of node N1. The settlement predictions of nodes N3 and N4 are approximately 8 and 4 cm, respectively. When nodes N1, N2, and N3 move horizontally toward the ocean (approximately 20 cm for nodes N1 and N2 and 8 cm for node N3), node N4 deforms to the right side of the model by around 2 cm. It is noticed that the displacements of soil behind the steel pipe piles (represented by node N4) are reduced compared to those of the soil in front of the row of steel piles (represented by nodes N1, N2, and N3).   Figure 22 illustrates the maximum shear strain contours at the end of the earthquake. Two shear bands are localized near the quay wall and the far field of the model. Near the quay wall, the shear band would be extended to the top of the quay wall due to the significant displacements of the caisson. Due to the free-field boundary conditions, the right side of the model was not constrained nor tied with a rigid body to fix the horizontal displacement. Therefore, many significant shear bands are localized in this region. Figure 23 depicts the shear stress-strain response during the earthquake of the observed elements in Figure 15. When maximum shear strain produced by the Mohr-Coulomb model is approximately 0.03% in element E4, the significantly higher ones generated by the UBCSAND2 model in elements E1, E2, and E3 are 0.35%, 0.2%, and 0.5%, respectively. Additionally, the effect of initial static shear stress (i.e., α > 0) can be addressed in the stress-strain responses of elements E1, E2, and E3 as shown in Figure 23a-c. This initial shear stress produced a constant bias in the cyclic shear load of the soil elements assigned by the UBCSAND2 model. In the regions behind the steel pipe piles represented by element E3, in which the α value is approximately zero, the symmetrical response of shear stress is in a range of (−30, 30) kPa. Whereas, the non-symmetrically mobilized shear stress of elements E1 with α ≈ 0.075 and E2 with α ≈ 0.05, are in the range of (−15, 25) and (−10, 15) kPa, respectively.  Figure 22 illustrates the maximum shear strain contours at the end of the earthquake. Two shear bands are localized near the quay wall and the far field of the model. Near the quay wall, the shear band would be extended to the top of the quay wall due to the significant displacements of the caisson. Due to the free-field boundary conditions, the right side of the model was not constrained nor tied with a rigid body to fix the horizontal displacement. Therefore, many significant shear bands are localized in this region.  Figure 22 illustrates the maximum shear strain contours at the end of the earthquake. Two shear bands are localized near the quay wall and the far field of the model. Near the quay wall, the shear band would be extended to the top of the quay wall due to the significant displacements of the caisson. Due to the free-field boundary conditions, the right side of the model was not constrained nor tied with a rigid body to fix the horizontal displacement. Therefore, many significant shear bands are localized in this region. Figure 23 depicts the shear stress-strain response during the earthquake of the observed elements in Figure 15. When maximum shear strain produced by the Mohr-Coulomb model is approximately 0.03% in element E4, the significantly higher ones generated by the UBCSAND2 model in elements E1, E2, and E3 are 0.35%, 0.2%, and 0.5%, respectively. Additionally, the effect of initial static shear stress (i.e., α > 0) can be addressed in the stress-strain responses of elements E1, E2, and E3 as shown in Figure 23a-c. This initial shear stress produced a constant bias in the cyclic shear load of the soil elements assigned by the UBCSAND2 model. In the regions behind the steel pipe piles represented by element E3, in which the α value is approximately zero, the symmetrical response of shear stress is in a range of (−30, 30) kPa. Whereas, the non-symmetrically mobilized shear stress of elements E1 with α ≈ 0.075 and E2 with α ≈ 0.05, are in the range of (−15, 25) and (−10, 15) kPa, respectively.    Figure 15. When maximum shear strain produced by the Mohr-Coulomb model is approximately 0.03% in element E4, the significantly higher ones generated by the UBCSAND2 model in elements E1, E2, and E3 are 0.35%, 0.2%, and 0.5%, respectively. Additionally, the effect of initial static shear stress (i.e., α > 0) can be addressed in the stress-strain responses of elements E1, E2, and E3 as shown in Figure 23a-c. This initial shear stress produced a constant bias in the cyclic shear load of the soil elements assigned by the UBCSAND2 model. In the regions behind the steel pipe piles represented by element E3, in which the α value is approximately zero, the symmetrical response of shear stress is in a range of (−30, 30) kPa. Whereas, the non-symmetrically mobilized shear stress of elements E1 with α ≈ 0.075 and E2 with α ≈ 0.05, are in the range of (−15, 25) and (−10, 15) kPa, respectively.  Figure 24 shows the comparison between the input ground motions and transmission motions in the analysis. The observed nodes were set in Figure 15. In this study, the recorded accelerationtime history of the Pohang earthquake was extracted to 40 s, filtered to 20 Hz, and baseline drifted. Then, a transformation of the corrected velocity motion to shear stress motion was made to be applicable to the quiet boundary condition in FLAC, as in Equation (19). Therefore, it is essential to check the wave transmission during shaking. At the bottom of the observed nodes, there are slight differences in the acceleration-time histories comparing the input ground motion. It can be said that the transformation from velocity history to shear stress history was successfully implemented. In dynamic analysis, it is important to consider the effect of velocity doubling at the model base and at the top of the model. If elastic materials are defined for the model and run in an un-damped analysis, the velocity-doubling effect would occur in the free surface [30]. In other words, the velocity at the top is double the velocity at the base for un-damped elastic materials. In this study, Rayleigh damping was used for the dynamic analysis, in which there is a small amplification of the velocity wave at the free surface compared to the base, indicating that the base is within the range of velocity. The trend of this effect is depicted in Figure 24b-d. It is interesting that the velocities at the top of the model for nodes N3 and N4 are approximately 25% higher than those at the base, when only 5% of the difference is presented for node N2. This may be due to the very high stiffness of the gravel preventing wave transmission from the base to the top surface of node N2.  Figure 24 shows the comparison between the input ground motions and transmission motions in the analysis. The observed nodes were set in Figure 15. In this study, the recorded acceleration-time history of the Pohang earthquake was extracted to 40 s, filtered to 20 Hz, and baseline drifted. Then, a transformation of the corrected velocity motion to shear stress motion was made to be applicable to the quiet boundary condition in FLAC, as in Equation (19). Therefore, it is essential to check the wave transmission during shaking. At the bottom of the observed nodes, there are slight differences in the acceleration-time histories comparing the input ground motion. It can be said that the transformation from velocity history to shear stress history was successfully implemented. In dynamic analysis, it is important to consider the effect of velocity doubling at the model base and at the top of the model. If elastic materials are defined for the model and run in an un-damped analysis, the velocity-doubling effect would occur in the free surface [30]. In other words, the velocity at the top is double the velocity at the base for un-damped elastic materials. In this study, Rayleigh damping was used for the dynamic analysis, in which there is a small amplification of the velocity wave at the free surface compared to the base, indicating that the base is within the range of velocity. The trend of this effect is depicted in Figure 24b-d. It is interesting that the velocities at the top of the model for nodes N3 and N4 are approximately 25% higher than those at the base, when only 5% of the difference is presented for node N2. This may be due to the very high stiffness of the gravel preventing wave transmission from the base to the top surface of node N2.

Conclusions
Since the widespread damages on geotechnical structures caused by liquefaction during the Pohang earthquake, the understanding of the mechanism and more accurate predictions of liquefaction-induced displacements using state-of-the-art numerical tools have played crucial roles in the design stages for practical engineering applications. In this study, a robustly practical model called UBCSAND2, in which two-mobilized planes (i.e., the maximum shear plane and horizontal plane) simultaneously participate in plastic strain generation, addresses the conventional effects in the element CDSS test: rotation of principal stress and plastic unloading.
The UBCSAND2 model was calibrated by a series of CDSS tests. There is a good agreement between the experimental and estimated results either in variations of CSR-Nliq curves or in terms of stress path, shear stress-shear strain, and pore pressure generation under different density states (Dr), confining stress (σ´v), and initial static shear stress (α) conditions. The generic parameters of the model were chosen to assess the liquefaction that occurred at Youngil-man port, Pohang, Korea.
The quay wall in the region of Pohang in South Korea was numerically simulated to demonstrate capabilities of the UBCSAND2 model in predicting the mechanical behavior of liquefied soils. According to the numerical results, liquefaction mostly occurred in silty sand layers when the excess pore water pressure generated almost 85−95% of initial vertical effective stress. The displacements at the top of the quay wall from liquefaction phenomenon were approximately 20 and 3 cm, respectively, when the settlements of reclaimed area behind the quay wall varied nearly from 4 to 8 cm. These numerical results obtained from the analysis are extremely similar to that observed by solely visiting the sites after the Pohang earthquake. The mechanism of excess pore water pressure generation, shear band expansion, shear stress-strain response, and transmission waves obtained from the numerical analysis have provided a full understanding of the liquefaction phenomenon that occurred in Youngil-man port.

Conclusions
Since the widespread damages on geotechnical structures caused by liquefaction during the Pohang earthquake, the understanding of the mechanism and more accurate predictions of liquefaction-induced displacements using state-of-the-art numerical tools have played crucial roles in the design stages for practical engineering applications. In this study, a robustly practical model called UBCSAND2, in which two-mobilized planes (i.e., the maximum shear plane and horizontal plane) simultaneously participate in plastic strain generation, addresses the conventional effects in the element CDSS test: rotation of principal stress and plastic unloading.
The UBCSAND2 model was calibrated by a series of CDSS tests. There is a good agreement between the experimental and estimated results either in variations of CSR-N liq curves or in terms of stress path, shear stress-shear strain, and pore pressure generation under different density states (D r ), confining stress (σ v ), and initial static shear stress (α) conditions. The generic parameters of the model were chosen to assess the liquefaction that occurred at Youngil-man port, Pohang, Korea.
The quay wall in the region of Pohang in South Korea was numerically simulated to demonstrate capabilities of the UBCSAND2 model in predicting the mechanical behavior of liquefied soils. According to the numerical results, liquefaction mostly occurred in silty sand layers when the excess pore water pressure generated almost 85−95% of initial vertical effective stress. The displacements at the top of the quay wall from liquefaction phenomenon were approximately 20 and 3 cm, respectively, when the settlements of reclaimed area behind the quay wall varied nearly from 4 to 8 cm. These numerical results obtained from the analysis are extremely similar to that observed by solely visiting the sites after the Pohang earthquake. The mechanism of excess pore water pressure generation, shear band expansion, shear stress-strain response, and transmission waves obtained from the numerical analysis have provided a full understanding of the liquefaction phenomenon that occurred in Youngil-man port.