Seismic Liquefaction of Saturated Calcareous Sands: Dynamic Centrifuge Test and Numerical Simulation

: Both a dynamic centrifuge test and dynamic ﬁnite element analysis were carried out to assess the seismic liquefaction risk of a saturated-calcareous-sand site in a port project in Timor-Leste. Taking the in situ calcareous sands as the model material, two groups of horizontal free-ﬁeld model tests for medium dense and dense saturated calcareous sands were completed based on the two-stage scaling law theory. Three natural earthquake records with varied peak accelerations were adopted as input motions. The experimental results indicate that the shallower the depth, the lower the relative density, the longer the seismic duration, the larger the peak acceleration, and the more susceptible the saturated site to liquefaction. The sands on site at a depth of ﬁve to ten meters is highly risky for liquefaction with the excess pore pressure ratio reaching up to about 1.0 under the seismic peak acceleration of 0.3 g. The risk of liquefaction for site sands is rather small under the seismic peak acceleration of 0.1 g. The study reveals the characteristics of the pore pressure development in sites of varied relative densities under different seismic loadings, which provides a scientiﬁc basis for the liquefaction risk assessment of the engineering site.


Introduction
Seismic-induced soil liquefaction can be destructive for infrastructures. History investigations have shown that seismic-induced soil liquefaction can cause large ground deformation or settlement, ground fissure, lateral spreading, etc. [1][2][3][4], as well as direct damage on structures including embankments, underground pipelines, and tunnels [5,6]. In recent years, catastrophes led by seismic-induced soil liquefaction took place worldwide, e.g., the 2008 Wenchuan earthquake, the 2011 Christchurch New Zealand earthquake, the 2011 east Japan earthquake, the 2016 southern Taiwan earthquake and the 2018 Sulawesi earthquake in Indonesia [7][8][9][10][11]. Significant soil liquefaction was observed during these earthquakes, which brought about destructive aftermath. Statistical investigations find that the economic loss caused by liquefaction can make up 15-30% of the total loss during an earthquake [12], which should be taken as one of the major hazards for earthquake damage.
Calcareous sand is widely spread in islands and coastal areas around the world. With the expansion of human society, the engineering projects confronting calcareous sand sites are increasing, and the seismic liquefaction risk of calcareous sand ground cannot be ignored. Differing from the commonly known quartzite sand, the calcareous sand possesses unique mechanical properties due to the distinct mesoscopic characteristics such as the angular grain shape, low stiffness, high particle breakage susceptibility, and high internal 2 of 17 porosity [13,14]. While considerable progresses have been made on studies about sand liquefaction, the majority of the previous studies are conducted based on quartzite sand [15]. Disputes and blanks still remain on the mechanical properties and dynamic behavior of calcareous sand. According to laboratory test results, it is believed that calcareous sand has a higher liquefaction resistance than quartzite sand, i.e., it is generally more difficult for calcareous sand to liquefy in most of the cases [16,17]. However, the history investigations can still present cases of seismic-induced liquefaction in calcareous sand sites, such as the 1993 Guam earthquake, the 2006 Hawaii earthquake, and the 2010 Haiti earthquake [18][19][20]. Sand liquefaction had caused severe damage to infrastructure including embankments and ports during these earthquakes [21]. Therefore, the liquefaction mechanism and dynamic characteristics of calcareous sand still need to be further studied.
The conventional methods of studying the seismic response of sites include history investigation, in situ testing, laboratory testing, and numerical simulation. Among the various types of laboratory tests, the centrifuge test stands out with the advantage of generating Ng gravity via rotation and obtaining an identical or approximately identical stress field to the prototype in the model, which makes it capable to reproduce the deformation and failure mechanism of the prototype [22]. The centrifuge tests have been widely used in earthquake geotechnical engineering studies for more than three decades, providing a large amount of essential information through physical models and greatly contributing to the advancement of knowledge in this field [23]. The pioneering collaborative project VELACS (Verification of Liquefaction Analysis by Centrifuge Studies) carried out in the mid-1990s provided reliable experimental data for numerically modelling soil liquefaction, including nine centrifuge model tests of saturated sand materials [24][25][26][27]. Another notable international collaboration, LEAP-2017 (Liquefaction Experiments and Analysis Projects), initiated for validating constitutive and numerical modelling of soil liquefaction and its consequences using centrifuge tests, eventually built a high-quality database including twenty-four centrifuge model tests of liquefaction and lateral spreading [28][29][30]. The LIQ-UEFACT project within the EU Horizon 2020 program just ending in October 2019 carried out a total of thirty-seven dynamic centrifuge tests including benchmark tests of saturated sand under free field conditions [31,32]. Extensive research other than theses remarkable collaborative projects has been conducted to evaluate soil liquefaction through centrifuge tests and could not be listed exhaustively.
On the other hand, based on the well-documented experimental data, the state-ofthe-art constitutive and numerical modelling methods for soil liquefaction are validated and further improved by researchers. Li et al. [33] assessed the seismic response of a layered liquefiable soil profile with granular columns through both numerical and centrifuge modeling. The numerical simulations are performed with PDMY and Manzari-Dafalias model in Opensees. Comparisons of numerical simulations with experimental results indicate a good prediction ability for both models on the peak value and generation of excess pore pressures during shaking. Similarly, implementing FE analyses with PDMY model in Opensees, Tang et al. [34] simulated the centrifuge test of the saturated silty site improved with granular columns and analyzed the dynamic characteristics such as ground acceleration and pore pressure. The results suggested a good agreement between the simulations with the experimental measurement. As one of the main work packages in LEAP-2017, eleven sets of Type-B predictions were implemented with PM4Sand model, PDMY, Hypoplastic model, Tsinghua constitutive model, etc. Manzari et al. [28] presents the comparisons of the Type-B simulations with the results of the centrifuge test and concluded that although similar trends can be achieved by a number of numerical simulations, predicting all the key quantities remains a challenge.
Ports are of high political and economic value, and their destructive consequences can be unbearable. The Tibar Bay Port Project (TBPP), located on the Circum-Pacific Belt, has a high seismic fortification intensity according to Chinese standards. The design seismic acceleration (475-year return period) of the wharf area reaches up to 0.53 g on the basis of the seismic hazard assessment report for TBPP. The site is mainly underlain by a mixed Ports are of high political and economic value, and their destructive consequences can be unbearable. The Tibar Bay Port Project (TBPP), located on the Circum-Pacific Belt, has a high seismic fortification intensity according to Chinese standards. The design seismic acceleration (475-year return period) of the wharf area reaches up to 0.53 g on the basis of the seismic hazard assessment report for TBPP. The site is mainly underlain by a mixed layer of sand and gravel with a maximum thickness of 40 m, which is very loose to loose, and a locally medium dense state.
A high proportion of fine sands can be found in local areas, and the average percentage of fines content smaller than 0.075 mm is 19% as shown in Figure 1. The layer has an average natural density of 1.65 g/cm 3 . The SPT value ranges from 0 to 18, with a mean value being 7, and the design relative density is 50%. The design permeability is 6.4 × 10 −2 cm/s. In order to scientifically evaluate the liquefaction potential of the calcareous sand site, centrifuge tests are carried out in the centrifuge shaking table system of Tongji University. Saturated calcareous sand models are tested with different relative densities under excitations of sets of input ground motions with different peak accelerations. The development characteristics of excess pore pressure at different depths during the earthquake are investigated, and the factors affecting excess pore pressure are analyzed in the tests. Moreover, based on the Biot's dynamic consolidation theory, seismic response analyses of the calcareous sand site are conducted on OpenSees. The dynamic pore pressure responses are simulated and compared with that in centrifuge models. The seismic liquefaction potential of calcareous sand sites is comprehensively evaluated.

Test Facilities
The dynamic centrifuge tests were conducted using the TLJ-150 geotechnical centrifuge with matching electro-hydraulic shaking table (see Figure 2) at Tongji University, China. The major specifications of the experimental facilities are shown in Table 1. The facilities can stably generate a maximum rotational acceleration of 50 g with a payload of 300 kg during the dynamic model tests, and a maximum shaking acceleration of 20 g can be produced with models weighing up to 180 kg. A laminar shear model box, which is composed of 24 stacked rings of rectangular-section hollow aluminum square tubes, was used in the tests. The stacked rings are connected by guide rails and can slide freely in the horizontal shaking direction. The maximum relative displacement between each stacked ring is 6 mm. The internal dimensions of the model box are 500 mm × 400 mm × 550 mm (length × width × height). A 1 mm-thick high-durability membrane is set inside the box to prevent leakage from the gap between the stacked rings. In order to scientifically evaluate the liquefaction potential of the calcareous sand site, centrifuge tests are carried out in the centrifuge shaking table system of Tongji University. Saturated calcareous sand models are tested with different relative densities under excitations of sets of input ground motions with different peak accelerations. The development characteristics of excess pore pressure at different depths during the earthquake are investigated, and the factors affecting excess pore pressure are analyzed in the tests. Moreover, based on the Biot's dynamic consolidation theory, seismic response analyses of the calcareous sand site are conducted on OpenSees. The dynamic pore pressure responses are simulated and compared with that in centrifuge models. The seismic liquefaction potential of calcareous sand sites is comprehensively evaluated.

Test Facilities
The dynamic centrifuge tests were conducted using the TLJ-150 geotechnical centrifuge with matching electro-hydraulic shaking table (see Figure 2) at Tongji University, China. The major specifications of the experimental facilities are shown in Table 1. The facilities can stably generate a maximum rotational acceleration of 50 g with a payload of 300 kg during the dynamic model tests, and a maximum shaking acceleration of 20 g can be produced with models weighing up to 180 kg. A laminar shear model box, which is composed of 24 stacked rings of rectangular-section hollow aluminum square tubes, was used in the tests. The stacked rings are connected by guide rails and can slide freely in the horizontal shaking direction. The maximum relative displacement between each stacked ring is 6 mm. The internal dimensions of the model box are 500 × 400 × 550 mm (length × width × height). A 1 mm-thick high-durability membrane is set inside the box to prevent leakage from the gap between the stacked rings. Appl

Scaling Law and Physical Models
To reproduce the same overburden stress in the scaled model as that in the prototype, the conventional scaling law for dynamic centrifuge tests is commonly established using the Buckingham π theorem. However, the shortcoming of the conventional scaling law is apparent: the set of scaling factors is built on the predetermined gravitational acceleration scaling factor with a relatively fixed relationship. Consequently, the large size engineering prototype could be difficult to model when the experiments are faced with restrictions of facility capacity. To overcome the limitation, Iai et al. proposed the concept of two-stage scaling and developed the generalized scaling law (GSL) in dynamic centrifuge tests [35]. As shown in Figure 3, μ is the dimensional scaling factor at the first stage, and η is the gravitational acceleration scaling factor at the second stage. The two-stage scaling and GSL extend the limits of the existing facilities and allow the dynamic centrifuge tests with a dimensional scaling factor reaching up to 100-1000 [35]. Due to the stress dependency of soil material, separate scales for stress and strain are needed to secure the kinetic and kinematic similarity between the model and the

Scaling Law and Physical Models
To reproduce the same overburden stress in the scaled model as that in the prototype, the conventional scaling law for dynamic centrifuge tests is commonly established using the Buckingham π theorem. However, the shortcoming of the conventional scaling law is apparent: the set of scaling factors is built on the predetermined gravitational acceleration scaling factor with a relatively fixed relationship. Consequently, the large size engineering prototype could be difficult to model when the experiments are faced with restrictions of facility capacity. To overcome the limitation, Iai et al. proposed the concept of two-stage scaling and developed the generalized scaling law (GSL) in dynamic centrifuge tests [35]. As shown in Figure 3, µ is the dimensional scaling factor at the first stage, and η is the gravitational acceleration scaling factor at the second stage. The two-stage scaling and GSL extend the limits of the existing facilities and allow the dynamic centrifuge tests with a dimensional scaling factor reaching up to 100-1000 [35].

Scaling Law and Physical Models
To reproduce the same overburden stress in the scaled model as that in the prototype, the conventional scaling law for dynamic centrifuge tests is commonly established using the Buckingham π theorem. However, the shortcoming of the conventional scaling law is apparent: the set of scaling factors is built on the predetermined gravitational acceleration scaling factor with a relatively fixed relationship. Consequently, the large size engineering prototype could be difficult to model when the experiments are faced with restrictions of facility capacity. To overcome the limitation, Iai et al. proposed the concept of two-stage scaling and developed the generalized scaling law (GSL) in dynamic centrifuge tests [35]. As shown in Figure 3, μ is the dimensional scaling factor at the first stage, and η is the gravitational acceleration scaling factor at the second stage. The two-stage scaling and GSL extend the limits of the existing facilities and allow the dynamic centrifuge tests with a dimensional scaling factor reaching up to 100-1000 [35]. Due to the stress dependency of soil material, separate scales for stress and strain are needed to secure the kinetic and kinematic similarity between the model and the Due to the stress dependency of soil material, separate scales for stress and strain are needed to secure the kinetic and kinematic similarity between the model and the prototype. Iai and Sungano [36] developed three different types of GSL (Iai's Type I, II, III GSL) within the framework forementioned, each based on a different strain scaling factor determination method. Iai's type II GSL assumes that the tangential stiffness of sand is proportional to the square root of the confining pressure, and the relations between different scaling factors is presented in Equations (1)- (3): where µ σ , µ ρ , µ D , µ ε are the scaling factors of stress, density, tangent stiffness, and strain, respectively. The applicability of the two-stage scaling has been investigated and validated in previous research. Tobita et al. [37] and Zhou et al. [38] adopted Iai's type II GSL in dynamic centrifuge tests of saturated sand and found that the model behaviors well mimic the prototype. In order to model the thick liquefiable calcareous sand layer in the engineering site, considering the limited amount of model material, the dimensions of the model container, and the capacity of the centrifuge, the two-stage scaling with Iai's Type II GSL was adopted in the dynamic centrifuge tests, the scaling factors of which are shown in Table 2. During the experiments, rotational acceleration was set to be 40 g, and µ and η were 2 and 40, respectively. The length, width, and depth of the saturated level ground model are 500 × 400 × 250 mm, corresponding to the prototype of 40 × 32 × 20 m. Additionally, two free-field models, medium dense (M1) and dense (M2), are tested.

Test Materials and Model Preparation
Horizontal free-field models of saturated calcareous sand were tested. The sand material is taken originally from the project site, and its grain size distribution is adjusted to align with that in Figure 1. The primary physical properties are shown in Table 3. It should be noted that the scaling law requires a lower permeability in the model than that in the prototype to reflect the actual accumulation and dissipation process of EXPP in the model. As indicated in Table 2, the permeability of the model is 1/µ 0.75 η of the prototype. It is realized in this study by using a carboxy methylcellulose (CMC) solution as a pore fluid substitute, which is transparent and has a close-to-water density. The CMC solution is obtained by mixing a predetermined amount of CMC powder with water to reach the desired viscosity. Therefore, the CMC solution with a viscosity of 67.3cst was adopted in the tests. The dry sand models were prepared using the compaction method by 4 layers. The mass and height of each layer is determined by the relative density, which is calculated through the maximum and minimum dry density. The dry sand model is then saturated in a vacuum container by flooding the fluid from the bottom up. In order to achieve a high saturation degree, the flow rate of CMC solution is controlled to be as slow as possible. After saturation, the height of the ground model was measured at 9 different positions on the surface, and the relative density of the model was obtained by calculating the mean value. The relative densities of the two models (M1 and M2) are 53.1% and 67.7%, respectively.

Instrumentation Layout
The study mainly focuses on the shaking acceleration of the sand, excess pore pressure, and ground settlement. Therefore, the instrumentations used include accelerometers, pore pressure transducers, and laser displacement transducers. As shown in Figure 4, the instrumentations are arranged by the following principles: (1) monitoring of the shaking table movement: accelerometers are placed on the shaking table surface to monitor the output motion of the table; (2) monitoring of the acceleration response of the sand: A row of single-axis accelerometers (A1-A5) are equally spaced on the left side of the model's central axis, and placed as close to the center as possible. (3) monitoring of pore pressure response in the sand: a row of pore pressure transducers (P1-P4) is equally spaced on the right side of model's central axis, and P1 , P3 are placed at the same depth as P1, P3, respectively, for verification; (4) monitoring of ground settlement: a laser displacement transducer D1 is arranged in the center of the model to monitor the settlement of the model. The dry sand models were prepared using the compaction method by 4 layers. The mass and height of each layer is determined by the relative density, which is calculated through the maximum and minimum dry density. The dry sand model is then saturated in a vacuum container by flooding the fluid from the bottom up. In order to achieve a high saturation degree, the flow rate of CMC solution is controlled to be as slow as possible. After saturation, the height of the ground model was measured at 9 different positions on the surface, and the relative density of the model was obtained by calculating the mean value. The relative densities of the two models (M1 and M2) are 53.1% and 67.7%, respectively.

Instrumentation Layout
The study mainly focuses on the shaking acceleration of the sand, excess pore pressure, and ground settlement. Therefore, the instrumentations used include accelerometers, pore pressure transducers, and laser displacement transducers. As shown in Figure  4, the instrumentations are arranged by the following principles: (1)

Input Ground Motions and Test Procedures
Three different earthquakes (E1, E2, and E3) defined in the seismic hazard assessment report for TBPP were selected as the input ground motions. These natural earthquake records have been adjusted to reflect the geological characteristics of the engineering site. The amplitude of acceleration for each input ground motion was linearly scaled to three levels of 0.1 g, 0.2 g, and 0.3 g, and the seismic loadings were input from the bottom of the model box along the lateral direction by a specified sequence as shown in Table 4. The white noise is applied to obtain any changes on the site's frequency spectrum

Input Ground Motions and Test Procedures
Three different earthquakes (E1, E2, and E3) defined in the seismic hazard assessment report for TBPP were selected as the input ground motions. These natural earthquake records have been adjusted to reflect the geological characteristics of the engineering site. The amplitude of acceleration for each input ground motion was linearly scaled to three levels of 0.1 g, 0.2 g, and 0.3 g, and the seismic loadings were input from the bottom of the model box along the lateral direction by a specified sequence as shown in Table 4. The white noise is applied to obtain any changes on the site's frequency spectrum characteristics in the course of being subject to seismic loadings with increasing peak acceleration amplitudes. An adequately long interval between every single seismic loading was needed for a full dissipation of the excess pore pressure, for which the subsequent seismic waves were not input until approximately 60 s (model time) after precedent excitation. The M1 and M2 models were applied with the same loading condition. Noted that the input motions of acceleration time histories must be scaled conforming to the scaling law. As indicated in Table 2, the time scaling factor is 1/67.3 (model/prototype, the same below), and the acceleration scaling factor is 40. Taking the actual input motions E1, E2, and E3 with a peak acceleration of 0. characteristics in the course of being subject to seismic loadings with increasing peak acceleration amplitudes. An adequately long interval between every single seismic loading was needed for a full dissipation of the excess pore pressure, for which the subsequent seismic waves were not input until approximately 60 s (model time) after precedent excitation. The M1 and M2 models were applied with the same loading condition. White noise 0.01 g Noted that the input motions of acceleration time histories must be scaled conforming to the scaling law. As indicated in Table 2, the time scaling factor is 1/67.3 (model/prototype, the same below), and the acceleration scaling factor is 40. Taking the actual input motions E1, E2, and E3 with a peak acceleration of 0.

Excess Pore Pressure Response along the Depth
This study mainly focuses on the seismic response of pore pressure to evaluate the liquefaction potential of the saturated calcareous sand site. All the experimental results in this paper have been transformed to the prototype values through the scaling law. The excess pore pressure ratio (EXPPR) is defined as the ratio of excess pore pressure to the effective overburden stress. The maximum values of EXPP and EXPPR along the depth of

Excess Pore Pressure Response along the Depth
This study mainly focuses on the seismic response of pore pressure to evaluate the liquefaction potential of the saturated calcareous sand site. All the experimental results in this paper have been transformed to the prototype values through the scaling law. The excess pore pressure ratio (EXPPR) is defined as the ratio of excess pore pressure to the effective overburden stress. The maximum values of EXPP and EXPPR along the depth of the M1 and M2 model subjected to different earthquakes is plotted in Figures 6 and 7. It is noted that the pore pressure time histories at a depth of 20 m for the M2 model was not completely recorded when subjected to the seismic loading with a peak acceleration of 0.3 g due to the damage. Moreover, the maximum values of EXPPR at the depth of 5 m and 10 m are just about 0.9 and 0.6, respectively, for M2 in the loading cases of E2−3 and E3−3. The lower EXPPR implies soils with a stronger liquefaction resistance than that in M1 and could be explained by the higher relative density of M2.

Excess Pore Pressure Time Histories during Strong Earthquakes
Further analysis is made on the EXPPR time histories during strong earthquakes. And the time histories of EXPPR at depths of 5 m, 10 m, and 15 m during the earthquakes Moreover, the maximum values of EXPPR at the depth of 5 m and 10 m are just about 0.9 and 0.6, respectively, for M2 in the loading cases of E2−3 and E3−3. The lower EXPPR implies soils with a stronger liquefaction resistance than that in M1 and could be explained by the higher relative density of M2.

Excess Pore Pressure Time Histories during Strong Earthquakes
Further analysis is made on the EXPPR time histories during strong earthquakes. And the time histories of EXPPR at depths of 5 m, 10 m, and 15 m during the earthquakes As shown in Figure 6, the seismic responses of pore pressure for M1 and M2 models behave quite differently. When subjected to different earthquakes, the maximum values of EXPP at different depths for both M1 and M2 models rise with the increase of the peak acceleration, and the increments get smaller in the deeper soils. The maximum values of EXPP resulted from earthquakes E2 and E3 are significantly higher than that of E1, and there is trivial difference between responses of E2 and E3. Moreover, compared with M1, the maximum values of EXPP for M2 are smaller at every depth.
Similar phenomena can be more explicitly observed in the maximum values of EXPPR along the depth. As shown in Figure 7, on the excitation of E2 and E3 with a peak acceleration of 0.3 g (E2−3 and E3−3), the maximum values of EXPPR at the depths of 5 m and 10 m both reached up to around 1.0 for M1, indicating that the soils have completely liquefied; the E1 with the same peak acceleration (E1−3) can only lead to peak values of 0.8 and 0.6 for EXPPR at the same depth of soils in M1; and the peak value of EXPPR at the depth of 20 m in M1 is only about 0.4 in the cases of E2−3 and E3−3. Generally, for different input motions, E2 and E3 cause a higher maximum value of pore pressure with the same peak acceleration; and when subjected to the same seismic loading, the EXPPR in the shallower soils present higher maximum values. The stronger liquefaction resistance of deeper soils might be relevant to the higher confining pressure, which brings about higher dynamic strength of soils and more significant contraction and densification behavior.
Moreover, the maximum values of EXPPR at the depth of 5 m and 10 m are just about 0.9 and 0.6, respectively, for M2 in the loading cases of E2−3 and E3−3. The lower EXPPR implies soils with a stronger liquefaction resistance than that in M1 and could be explained by the higher relative density of M2.

Excess Pore Pressure Time Histories during Strong Earthquakes
Further analysis is made on the EXPPR time histories during strong earthquakes. And the time histories of EXPPR at depths of 5 m, 10 m, and 15 m during the earthquakes of E1−3, E2−3, and E3−3 is plotted for both M1 and M2, which is shown in Figure 8.
It can be seen that for the same seismic loading and the same model, the peaks of the EXPPR of the soils decrease conspicuously with the increase of depth, and meanwhile the dissipation of pore pressure slightly accelerates; still for the same seismic loading, the peaks of EXPPR in M1 are significantly higher than that in M2. However, the dissipation rates of M1 and M2 are similar, with M2 being slightly faster in most of the cases. Moreover, it is worth noticing that during the earthquakes, the time histories of EXPPR in M2 exhibits more obvious "spikes", which indicates that the high relative density calcareous sand of M2 is more easily able to dilate under dynamic loading. Among the earthquakes of E1, E2, and E3, E3 leads to the fastest accumulation of EXPP in the soils, while E1 leads to the slowest.
Apart from the aforementioned characteristic, it should also be noted that under the same seismic loading, the time histories of EXPPR at a depth of 15 m display significantly smaller "spikes" for both M1 and M2. The phenomenon might be explained by the fact that the soils under high confining pressure have a relatively stronger tendency toward contraction rather than dilation during the dynamic loading, which agrees with the discussion in last section.
peaks of EXPPR in M1 are significantly higher than that in M2. However, the dissipation rates of M1 and M2 are similar, with M2 being slightly faster in most of the cases. Moreover, it is worth noticing that during the earthquakes, the time histories of EXPPR in M2 exhibits more obvious "spikes", which indicates that the high relative density calcareous sand of M2 is more easily able to dilate under dynamic loading. Among the earthquakes of E1, E2, and E3, E3 leads to the fastest accumulation of EXPP in the soils, while E1 leads to the slowest. Apart from the aforementioned characteristic, it should also be noted that under the same seismic loading, the time histories of EXPPR at a depth of 15 m display significantly smaller "spikes" for both M1 and M2. The phenomenon might be explained by the fact that the soils under high confining pressure have a relatively stronger tendency toward contraction rather than dilation during the dynamic loading, which agrees with the discussion in last section.

Constitutive Model and Parameters
The centrifuge tests could be inevitably restricted by realistic factors such as facility capacity, testing materials, and experiment costs, which pose difficulties in physically modelling the complex boundary value problems in engineering. With advantages in this aspect, the dynamic numerical analysis has become an important alternative to study the seismic response of saturated sites. Additionally, the FE analysis based on Biot's dynamic consolidation equation, which fully couples the pore pressure with the dynamic response of the soil skeleton, is one of the most promising and renowned numerical approaches [39], within the framework of which the constitutive model of soil plays a critical role and determines its prediction capability. The multi-yield surface elastoplastic constitutive, PDMY02, was applied in this study to simulate the liquefiable calcareous sand material [40,41]. Additionally, the OpenSees, an open-source computation platform, was used to carry out the fully coupled FE analyses on the seismic response of saturated ground. The prediction capability of the numerical approach was also evaluated.
PDMY can well describe the cyclic mobility and irreversible shear strain accumulation behavior of liquefiable sand under seismic excitation, which has been verified in many research applications [33,34,[42][43][44][45]. It is developed based on the multiple yield surface theory proposed by Prevost [46], and a combination of associative and non-associative flow rules is adopted. A total of 17 input parameters are needed for the model. Gingery et al. [47] calibrated the parameters of sand with different densities, the method of which is referred to in this paper for calibrating the sand materials, combined with the insitu measurement. The parameters for constitutive model are listed in Table 5.

Constitutive Model and Parameters
The centrifuge tests could be inevitably restricted by realistic factors such as facility capacity, testing materials, and experiment costs, which pose difficulties in physically modelling the complex boundary value problems in engineering. With advantages in this aspect, the dynamic numerical analysis has become an important alternative to study the seismic response of saturated sites. Additionally, the FE analysis based on Biot's dynamic consolidation equation, which fully couples the pore pressure with the dynamic response of the soil skeleton, is one of the most promising and renowned numerical approaches [39], within the framework of which the constitutive model of soil plays a critical role and determines its prediction capability. The multi-yield surface elastoplastic constitutive, PDMY02, was applied in this study to simulate the liquefiable calcareous sand material [40,41]. Additionally, the OpenSees, an open-source computation platform, was used to carry out the fully coupled FE analyses on the seismic response of saturated ground. The prediction capability of the numerical approach was also evaluated.
PDMY can well describe the cyclic mobility and irreversible shear strain accumulation behavior of liquefiable sand under seismic excitation, which has been verified in many research applications [33,34,[42][43][44][45]. It is developed based on the multiple yield surface theory proposed by Prevost [46], and a combination of associative and non-associative flow rules is adopted. A total of 17 input parameters are needed for the model. Gingery et al. [47] calibrated the parameters of sand with different densities, the method of which is referred to in this paper for calibrating the sand materials, combined with the in-situ measurement. The parameters for constitutive model are listed in Table 5.   − 2ϑ)). The Poisson's ratio is set to be 0.499 and 0, respectively, for the stress field initialization and dynamic computation.

Finite Element Model and Computation Procedure
The 2D FE model of the saturated calcareous sand with relative density of 50% was established and analyzed. As shown in Figure 9, the dimensions of the FE model are 20 m × 40 m (height and width), the height of which is consistent with that of the centrifuge model. The sand material is modeled by QuadUP elements [48] following the u-p dynamic consolidation equation, and the thickness of elements is 1 m. According to the laboratory test results, the permeability k sand of sand materials is set to be 6.4 × 10 −2 cm/s. The boundary conditions are set as follows: (1) the displacement of nodes in left and right boundaries sharing the same elevation are forced to be equal in x-direction and y-direction; (2) The nodes at the base are fixed, and the input seismic loading is imposed on the base along the x-direction; (3) The pore pressure on the soil surface is fixed to be zero, and other boundaries are impervious. The earthquake E2 in the centrifuge tests are selected as the input seismic loading, and the peak acceleration is linearly scaled to 0.1 g, 0.2 g, and 0.3 g.  1 2 )). The Poisson's ratio is set to be 0.499 and 0, respectively, for the stress field initialization and dynamic computation.

Finite Element Model and Computation Procedure
The 2D FE model of the saturated calcareous sand with relative density of 50% was established and analyzed. As shown in Figure 9, the dimensions of the FE model are 20 m × 40 m (height and width), the height of which is consistent with that of the centrifuge model. The sand material is modeled by QuadUP elements [48] following the -dynamic consolidation equation, and the thickness of elements is 1 m. According to the laboratory test results, the permeability of sand materials is set to be 6.4 × 10 −2 cm/s. The boundary conditions are set as follows: (1) the displacement of nodes in left and right boundaries sharing the same elevation are forced to be equal in -direction and -direction; (2) The nodes at the base are fixed, and the input seismic loading is imposed on the base along the -direction; (3) The pore pressure on the soil surface is fixed to be zero, and other boundaries are impervious. The earthquake E2 in the centrifuge tests are selected as the input seismic loading, and the peak acceleration is linearly scaled to 0.1 g, 0.2 g, and 0.3 g. Figure 9. Diagram of FE model for level ground of saturated calcareous sand. Figure 9. Diagram of FE model for level ground of saturated calcareous sand.
Uniform linear elastic mechanical properties were applied to all the materials in the model to establish the initial static stress condition as well as to avoid the stress concentration or arching effects. After that, the corresponding input parameters of sand materials were restored, and elastoplastic analysis was performed to obtain the initial stress field and pore pressure field. The ground motion was imposed on the base of the model by means of acceleration time history, namely the vibration method. During the computation, Rayleigh damping of 0.5% was applied to control high-frequency noise and ensure the damping energy dissipation of the soils under small strain condition.

Simulation Results
Observation points A, B, and C were placed sequentially at a depth of 5 m, 10 m, and 15 m along the central axis of the model, as indicated in Figure 9. The seismic response of EXPPR at the observation points during earthquakes of different peak accelerations were plotted and compared with the measured results in the centrifuge model (see Figure 10).
Uniform linear elastic mechanical properties were applied to all the materials in the model to establish the initial static stress condition as well as to avoid the stress concentration or arching effects. After that, the corresponding input parameters of sand materials were restored, and elastoplastic analysis was performed to obtain the initial stress field and pore pressure field. The ground motion was imposed on the base of the model by means of acceleration time history, namely the vibration method. During the computation, Rayleigh damping of 0.5% was applied to control high-frequency noise and ensure the damping energy dissipation of the soils under small strain condition.

Simulation Results
Observation points A, B, and C were placed sequentially at a depth of 5 m, 10 m, and 15 m along the central axis of the model, as indicated in Figure 9. The seismic response of EXPPR at the observation points during earthquakes of different peak accelerations were plotted and compared with the measured results in the centrifuge model (see Figure 10).  With regard to the numerical results, it can be seen that with the increase of the peak acceleration of the input ground motion, the EXPPRs at all the depths reach up to a higher maximum value with faster accumulation and slower dissipation. Under the earthquake with peak acceleration of 0.3 g, the EXPPRs at all the observation points reach 1.0 after about 9 s, even for point C located at 15 m depth, implying a high risk of liquefaction for the simulated site. Moreover, the EXPPRs remain in a high platform and do not dissipate during the loading of the earthquake. The earthquakes with peak acceleration of 0.1 g and 0.2 g both result in EXPPRs at different depths lower than 0.6, indicating that the liquefaction potential of the site is relatively low.
Compared with the experimental results, the numerical results generally present a faster process in both accumulation and dissipation but relatively consistent peak values of the EXPPRs. The only noticeable discrepancy occurs at the depth of 15 m under PGA of 0.3 g. The numerical results largely overestimated the EXPPR, which may be attributed to the overburden stress effect that is not well incorporated in the constitutive model. Moreover, it is obvious that the numerical results hardly have the "spikes", which is distinct in experiments. Generally, both the numerical and experimental results imply that the saturated calcareous sand site with relative density of 50% will have a high liquefaction potential for the sand at shallow depth when subjected to the earthquakes having peak acceleration of 0.3 g or even higher. In addition, the dynamic FE analysis framework adopting -formulation and PDMY02 was proved to effectively predict the peak values of EXPPR during earthquake.

Conclusions
In this paper, two-stage scaling law theory was adopted to carry out dynamic centrifuge tests on two calcareous sand models with different relative densities (M1 and M2) under three different earthquakes by means of stepwise loading. Moreover, dynamic FE analyses were conducted on a medium dense site subjected to strong earthquakes. The dynamic responses of the horizontal saturated free field during earthquakes were analyzed with a focus on the characteristics of pore pressure development. The liquefaction potential of the engineering site is comprehensively evaluated, and three major conclusions are drawn as follows: With regard to the numerical results, it can be seen that with the increase of the peak acceleration of the input ground motion, the EXPPRs at all the depths reach up to a higher maximum value with faster accumulation and slower dissipation. Under the earthquake with peak acceleration of 0.3 g, the EXPPRs at all the observation points reach 1.0 after about 9 s, even for point C located at 15 m depth, implying a high risk of liquefaction for the simulated site. Moreover, the EXPPRs remain in a high platform and do not dissipate during the loading of the earthquake. The earthquakes with peak acceleration of 0.1 g and 0.2 g both result in EXPPRs at different depths lower than 0.6, indicating that the liquefaction potential of the site is relatively low.
Compared with the experimental results, the numerical results generally present a faster process in both accumulation and dissipation but relatively consistent peak values of the EXPPRs. The only noticeable discrepancy occurs at the depth of 15 m under PGA of 0.3 g. The numerical results largely overestimated the EXPPR, which may be attributed to the overburden stress effect that is not well incorporated in the constitutive model. Moreover, it is obvious that the numerical results hardly have the "spikes", which is distinct in experiments. Generally, both the numerical and experimental results imply that the saturated calcareous sand site with relative density of 50% will have a high liquefaction potential for the sand at shallow depth when subjected to the earthquakes having peak acceleration of 0.3 g or even higher. In addition, the dynamic FE analysis framework adopting u-p formulation and PDMY02 was proved to effectively predict the peak values of EXPPR during earthquake.

Conclusions
In this paper, two-stage scaling law theory was adopted to carry out dynamic centrifuge tests on two calcareous sand models with different relative densities (M1 and M2) under three different earthquakes by means of stepwise loading. Moreover, dynamic FE analyses were conducted on a medium dense site subjected to strong earthquakes. The dynamic responses of the horizontal saturated free field during earthquakes were analyzed with a focus on the characteristics of pore pressure development. The liquefaction potential of the engineering site is comprehensively evaluated, and three major conclusions are drawn as follows: (1) The maximum value of EXPPR for the saturated calcareous sand ground will decrease with the increase of soil's depth, while it will increase with the increase of the input ground motion's peak acceleration and duration. The improvement of the relative density of can effectively strengthen the liquefaction resistance of the ground, reduce the peak pore pressure, and accelerate the dissipation process. (2) When subjected to earthquakes with 0.3 g peak acceleration, the saturated calcareous sand ground model witnessed the peak values of EXPPR reaching up to around 1.0 at the depth of 5-10 m, indicating a high potential of liquefaction under the same circumstances. (3) Based on Biot's dynamic consolidation theory, the dynamic FE analysis using the PDMY constitutive model can well simulate the seismic response of the EXPP at the initial stage of pore pressure accumulation. However, generally, the pore pressure could build up and dissipate relatively faster in the numerical results. Institutional Review Board Statement: Not applicable.

Informed Consent Statement: Not applicable.
Data Availability Statement: Not applicable.