Liquefaction analysis of marine diversion dike foundation and gravel pile reinforcement treatment based on the PL-Finn model

TL;DR

This study examines how gravel piles can prevent soil liquefaction in coastal nuclear power plant foundations during earthquakes.

Contribution

The study introduces a 3D numerical model using FLAC3D and the PL-Finn model to evaluate gravel pile reinforcement against liquefaction.

Findings

Unreinforced foundations show significant residual deformation from liquefaction.

Gravel pile reinforcement reduces excess pore water pressure and improves stability.

Gravel piles effectively mitigate liquefaction under design-level seismic intensity.

Abstract

This study investigates the liquefaction characteristics and failure mechanisms of sand foundations in the marine diversion dike of a coastal nuclear power plant, along with the anti-liquefaction performance of gravel pile-reinforced foundations. Using the three-dimensional finite difference software FLAC3D (Fast Lagrangian Analysis of Continua), we conducted numerical simulations of both unreinforced and reinforced liquefiable sand foundations. A three-dimensional model was developed, incorporating dynamic analysis with the PL-Finn constitutive model to simulate post-liquefaction large deformations. The efficacy of compacted gravel piles was evaluated through residual deformation, excess pore water pressure, and pore pressure ratio. Results demonstrate that unreinforced foundations exhibit systematic residual deformation due to liquefaction-induced sand flow, which is significantly reduced by gravel pile reinforcement. Both excess pore water pressure and pore pressure ratio decrease markedly after reinforcement. This substantially enhances the composite foundation’s structural stability and liquefaction resistance. These findings confirm that compacted gravel pile groups provide effective drainage and anti-liquefaction capacity, meeting site liquefaction mitigation requirements under design-level ground motion intensity, thereby offering a reference for similar engineering projects.

Figures

Diagram of diversion dike structure and foundation cross-section.

Three-dimensional finite element model of diversion dike.

Three-dimensional finite element model of the diversion dike after reinforcement.

Schematic diagram of gravel pile group arrangement (a) Plan view schematic (b) Pile layout schematic.

Seismic wave time-history curve (a) horizontal direction 1 (b) horizontal direction 2 (c) vertical direction.

The schematic diagram of free-field boundary.

Schematic diagram of monitoring points and units.

The analysis procedure.

Cloud chart of residual deformation distribution before reinforcement (a) Cloud chart of horizontal residual deformation distribution (b) Cloud chart of Vertical residual deformation distribution.

Cloud chart of residual deformation distribution after reinforcement (a) Cloud chart of horizontal residual deformation distribution (b) Cloud chart of Vertical residual deformation distribution.

Time-history curve of deformation at monitoring points before and after reinforcement (a) Horizontal displacement deformation curve of point a (b) Vertical displacement deformation curve of point a (c) Horizontal displacement deformation curve of point b (d) Vertical displacement deformation curve of point b.

Time-history curves of excess pore water pressure for different monitoring units: (a) Element A (b) Element B (c) Element C (d) Element D (e) Element E.

Time-history curves of excess pore pressure ratio for different monitoring units: (a) Element A; (b) Element B; (c) Element C; (d) Element D; (e) Element E.

Distribution of liquefied zones in the three-dimensional model: (a) Before reinforcement; (b) After reinforcement.

Equations

Full text

1 Introduction

Historical seismic damage records indicate that structural failures in diversion dikes are frequently caused by foundation sand liquefaction. Consequently, foundation treatment for diversion dike structures constitutes a critical aspect of earthquake and liquefaction resistance. For many years, researchers worldwide have focused on enhancing the liquefaction resistance of diversion dikes, with the goal of mitigating structural damage and losses during earthquakes. Numerous scholars have conducted extensive model testing and numerical analysis to investigate the seismic performance of these structures, yielding significant findings. Yang et al. [1] employed LS-DYNA software to perform dynamic elastic-plastic analysis, incorporating soil-structure interaction (SSI). They analyzed the seismic response of a nuclear power plant seawall and the influence of seismic wave parameters, providing theoretical support for its design and seismic performance. Liu et al. [2] studied the seismic response of the diversion dike head at a coastal nuclear power plant. They established a three-dimensional model using ANSYS and employed FLAC^3D^ coupled with the PL-Finn model for static and dynamic analysis. Zhao et al. [3] utilized the Hardin-Drnevich equivalent linear model to establish a two-dimensional model of a nuclear power plant diversion dike. Using a dynamic effective stress analysis method, they calculated the stability safety factor, residual deformation, and liquefaction area under SL1 and SL2 earthquake conditions.

Research on soil-structure dynamic interaction (SSI) originated from the study by Reissner et al. [4] on vibrations of rigid circular foundations on an elastic half-space surface. With advances in modern numerical methods and computer technology, as well as the successive construction of major engineering projects, SSI research under seismic loads has progressed significantly in both theory and practice, yielding substantial advances [5–9].

By reinforcing liquefiable soil layers with gravel piles to form a composite foundation, the bearing capacity is enhanced, settlement deformation reduced, and excess pore water pressure dissipation accelerated due to high permeability [10]. Gravel piles are extensively adopted in engineering practice owing to readily available materials, rapid construction, and cost-effectiveness [11]. Despite these applications, theoretical understanding—particularly regarding seismic liquefaction resistance—lags behind. Current findings primarily rely on field tests (e.g., in-situ and geotechnical investigations), yet complex engineering scenarios often cannot be fully resolved experimentally. Therefore, numerical simulations should complement experimental analyses where feasible, validating and supplementing field data to provide actionable guidance for on-site testing. To simulate liquefaction resistance of gravel pile-reinforced foundations under parametric variations, researchers have employed diverse three-dimensional finite-difference software. Jiang et al. [12] numerically analyzed deformation in expressway gravel pile composite foundations under varying conditions: pile construction processes, layered embankment/pavement filling, and traffic loads. Using FEMEPDYN, Li [13] proved that arranging gravel piles with long-short alternation (long piles on sides, short in center) enhances both liquefaction resistance and cost-efficiency in sand foundations. Huang [14] systematically reviewed prior studies and concluded that gravel piles provide comparable liquefaction resistance in silty and sandy foundations. Zhang [15] compared concrete piles versus gravel piles for reinforcing inclined liquefiable sites, finding concrete piles superior under identical embedment depths and pile diameters. Through numerical simulation, Shen [16] demonstrated that extending gravel pile lengths in saturated silt layers significantly reduces vertical/horizontal displacements of subway shield tunnels. For the Sihu Expressway’s silty sand foundation, Wang [17] implemented gridded gravel piles. Experimental tests and FLAC^3D^ analyses demonstrated their dual role: drainage/pressure relief for liquefaction mitigation and enhanced shear capacity, proving highly effective against seismic liquefaction in saturated silt. Gu [18] investigated gravel piles as anti-liquefaction measures for bridge foundations, detailing design methodologies and validation approaches. Niu et al. [19] developed a full soil-structure model in FLAC^3D^ to quantitatively assess soil-cement pile efficacy via residual deformation, excess pore pressure (ratio), and surface acceleration. Liu [20] simulated single versus group gravel piles reinforcing liquefiable sandy foundations via FEM, revealing superior performance of pile groups. Using FLAC^3D^’s Finn module, Li [21] confirmed that gravel piles significantly improve dynamic responses in both structures and soil layers. Qiu [22] proposed and validated an OpenSees-based evaluation framework for assessing liquefaction resistance in gravel pile-reinforced saturated sand foundations, corroborated by model tests. Through three-dimensional finite element simulations in OpenSeesPL, Asgari et al. [23] parametrically analyzed gravel piles and dowel pins for liquefaction mitigation. Key factors included area replacement ratio, soil/pile permeability, slope angle, pile diameter, superstructure mass, and seismic characteristics.

In this study, we establish dynamic analysis models of composite foundations for the saturated sand foundation of a marine diversion dike at a coastal nuclear power plant, comparing pre- and post-gravel pile construction conditions using FLAC^3D^ (Fast Lagrangian Analysis of Continua). Employing the PL-Finn constitutive model [24] for post-liquefaction large deformation, we conduct dynamic analyses to evaluate the anti-liquefaction efficacy of compacted gravel piles through three key indicators: residual deformation, excess pore water pressure, and pore pressure ratio. This approach enables comprehensive assessment of the anti-liquefaction properties, mechanisms, and reinforcement effects of gravel pile composite foundations, establishing a basis for optimizing gravel pile foundation reinforcement schemes.

2 Principles of calculation and analysis

2.1 Principle of FLAC3D dynamic analysis

FLAC^3D^ adopts a completely nonlinear dynamic response analysis method, which can accurately simulate plastic collapse and plastic flow characteristics of geomaterials. Besides, in the dynamic analysis module, the user-defined nonlinear constitutive model can be used to reflect the stress-strain hysteresis and deformation characteristics of soil. In addition, FLAC^3D^ can also adopt a fully coupled fluid-mechanics analysis model, in which the pore pressure response of the element affects the plastic volume deformation, and the plastic volume deformation of the element affects the pore water pressure. It can simulate the law of pore pressure growth and dissipation during the sand liquefaction process. The basic controlling formula of liquid-solid coupling in FLAC^3D^ dynamic calculation is described by the following fluid continuity formula:

$$[eqn]$$

In the above formula, $[eqn]$ is the impulse boundary, which is defined as the volume change per unit volume of fluid in porous media. $[eqn]$ represents body strain; the calculation size is related to the elastoplastic constitutive formula. M is the Biot modulus and $[eqn]$ is the Biot coefficient.

2.2 Liquefaction discrimination

Under dynamic action such as earthquakes, the saturated foundation will be changed from solid to liquid during liquefaction. When the viscous force of liquid is not considered, its shear strength will be changed to 0. In other words, when effective principal stress is 0, the saturated foundation will be liquefied. According to the principle of effective stress, the principal stress in three directions of the soil should be equal to the pore water pressure. Therefore, $[eqn]$ . $[eqn]$ is the principal stress in three directions of soil during liquefaction, kpa; $[eqn]$ is the pore water pressure during liquefaction, kpa.

The excess pore water pressure ratio $[eqn]$ can be expressed as:

$$[eqn]$$ $$[eqn]$$ $$[eqn]$$

$[eqn]$ is the average effective stress of the unit before the dynamic calculation, kpa and $[eqn]$ is the average effective stress of the unit during the dynamic calculation, kpa.

2.3 Fluid constitutive model after liquefaction

The constitutive model of sand flow characteristics after liquefaction is mainly established according to two different states of sand, namely zero effective stress state and non-zero effective stress state. The shear stress-shear strain rate relationship under zero effective stress state is as follows:

$$[eqn]$$

The shear stress-shear strain rate relationship under a non-zero effective stress state is:

$$[eqn]$$

In this formula, $[eqn]$ represents the generalized shear stress, kpa. $[eqn]$ represents the generalized strain rate; $[eqn]$ , $[eqn]$ are the fitting parameters, and $[eqn]$ is the excess pore water pressure ratio.

The Finn model takes deformation reversal into account and analyzes the problem in three dimensions, with at least six strain rate tensor components. These six strain rate tensors are defined as:

$$[eqn]$$ $$[eqn]$$

In this formula, e_ij_ is the strain increment tensor.

The extreme point of its trajectory is found in the strain space. (^o^) represents the previous point, (^oo^) represents the second point, and the previous unit vector n^o^ in the strain space is calculated as follows:

$$[eqn]$$ $$[eqn]$$

In this formula, v_i_ represents the strain increment on each component, and the subscript i represents 1–6 components.

New $[eqn]$ projection d onto the previous unit vector d is:

$$[eqn]$$

If d < 0, the strain is reversed. The absolute value of d is monitored during calculation. When the minimum calculation period is reached, if |d| reaches the maximum value |d|>=d max, then the inverse calculation is carried out, that is:

$$[eqn]$$

After shear strain $[eqn]$ and plastic volume strain increment $[eqn]$ are obtained, the plastic volumetric strain is updated:

$$[eqn]$$

According to $[eqn]$ , the shear strain increment can be modified:

$$[eqn]$$

In FLAC^3D^, the compressive strain increment has negative values and the plastic volume strain has positive values. Therefore, the average effective stress gradually decreases with the increase of plastic volume strain accumulation.

Based on the original liquefaction calculation model (the Finn model), the professional software of geotechnical engineering FLAC^3D^ has improved and added the post-liquefaction flow calculation module to form the PL-Finn model. It can more comprehensively consider the initial consolidation state and stress-strain relationship of soil before liquefaction, as well as the residual stress and residual deformation of soil after liquefaction. The PL-Finn model development procedure diagram is illustrated in S1 Fig.

2.4 Liquefaction mechanism of foundation elimination by gravel piles

The principle of vibration gravel piles to reinforce soft soil foundation is mainly to replace part of soft soil, form holes in sand, silty soils, cohesive soil, artificial fill, miry soil and other soil layers, and then backfill gravel and other coarse materials to form pile. Moreover, form a composite foundation with the original foundation soil to meet the bearing capacity required by the project, meet the settlement less than the limit value of the project, reduce deformation and compressibility and eliminate soil liquefaction. The principle of anti-liquefaction of gravel piles is mainly divided into three aspects: (1) compacting effect; (2) drainage effect; (3) pre-vibration effect. The bearing capacity of the gravel pile composite foundation is composed of the bearing capacity characteristic value f_spk_ of the gravel piles in the foundations and the bearing capacity characteristic value f_sk_ of soil between piles.

3 Engineering cases

3.1 Engineering situation

The marine engineering infrastructure for a nuclear power plant primarily comprises diversion dikes and revetments. This study focuses on a typical cross-section of the diversion dike. The site’s upper strata consist of Quaternary deposits (marine, alluvial, and residual), underlain by Mesozoic bedrock. From top to bottom, the foundation soils include: mucky soil, silty clay, medium-coarse sand, fully weathered granite, highly weathered granite, and moderately weathered granite. Based on local hydrogeological data, the design water level is set at −2.23 m below datum. Fig 1 illustrates the diversion dike structure and foundation cross-section.

Diagram of diversion dike structure and foundation cross-section.

3.2 Calculation models and parameters

Fig 2 shows the three-dimensional finite element model of the diversion dike. The model is 300m long, 100m wide and 55m high. It is composed of 134,242 units and 144,017 nodes. The unit size is 1.5-2.5m. In the analysis, the medium-coarse sand layer between silty clay and fully weathered rock is set as the liquefaction layer; the PL-Finn model is used for medium-coarse sand layer, and the Mohr-Coulomb model is used for other soil layers.

Three-dimensional finite element model of diversion dike.

Fig 3 is a three-dimensional finite element model of the diversion dike after reinforcement. The pile structural unit of FLAC^3D^ defines the compacted gravel pile group. Several scholars analyzed the influence range of pile diameter, pile length and pile spacing on pore water pressure, compaction degree and soil around piles by means of model tests and numerical simulation, and they proposed the reasonable pile spacing should be 2.5d ~ 3d [25]. Therefore, the foundation reinforcement adopts an equilateral triangle pile layout, in which the gravel pile diameter is 1.2 m, gravel pile spacing is 3.6m, and 450 piles in total. The length of the gravel pile is set at 25m to penetrate the mucky layer, silty clay layer and medium-coarse sand layer. A gravel cushion layer with 0.5m thickness is laid between the top of the compacted gravel pile and the diversion dike. Fig 4 shows the layout plan of the gravel pile.

Three-dimensional finite element model of the diversion dike after reinforcement.

Schematic diagram of gravel pile group arrangement (a) Plan view schematic (b) Pile layout schematic.

The mechanical parameters of rock-soil mass, fluid parameters and liquefaction parameters in the analysis are shown in [Tables 1–3](#pone.0330325.t001 pone.0330325.t002 pone.0330325.t003). The parameters are selected according to the Geotechnical Engineering Investigation Report and Test Report on Static and Dynamic Characteristics of Rock and Soil of the nuclear power plant. The soil mechanics parameters and contact surface parameters of the compacted gravel piles are shown in Table 4 and 5 below. The research direction of the contact surface is the “soft” contact surface. Nuclear power engineering adopts a conservative design, and it is necessary to consider the ultimate load-bearing condition when the contact surface is in a completely smooth state. At this time, the anti-slip ability of the structure is the weakest. Thus, specific friction angle parameters are assigned to the contact surface.

Table 1: Residual deformation parameters of soil.

Table 2: Soil mechanical parameters.

Table 3: Soil fluid and liquefaction parameters.

Table 4: Soil mechanical parameters of compacted gravel piles.

Table 5: Parameters of the contact surface.

The bulk modulus and shear modulus are required for calculation and analysis in FLAC^3D^, but the above parameters are not provided. Therefore, calculate according to the following formula:

$$[eqn]$$ $$[eqn]$$

In this formula, $[eqn]$ is the bulk modulus, $[eqn]$ is the shear modulus, $[eqn]$ is the Yang’s modulus, $[eqn]$ is the Poisson ratio.

When doing fluid calculations in FLAC^3D^, the permeability coefficient $[eqn]$ is one of the key parameters. Usually, the permeability coefficient $[eqn]$ of soil mechanics provided in the investigation report cannot be directly applied to FLAC^3D^. It needs to be converted by the formula in the calculation.

$$[eqn]$$

3.3 The setting of seismic dynamic loads and dynamic boundary conditions

According to the seismic safety assessment materials of a coastal nuclear power plant, the dynamic time-history analysis adopts the improved RG1.60 response spectrum of the United States, whose damping ratio is 5% and the duration is 25s. The corresponding time-history curve of seismic waves is shown in Fig 5. The design value of the horizontal peak acceleration of the bedrock is 0.30g, and the vertical peak acceleration is 0.20g. The free field boundary is set around the dynamic calculation model, and local damping is set to make the finite element mesh of the main structure coupled with the free field mesh, and the damping coefficient is 0.157. The boundary condition schematic diagram is shown in Fig 6.

Seismic wave time-history curve (a) horizontal direction 1 (b) horizontal direction 2 (c) vertical direction.

The schematic diagram of free-field boundary.

3.4 Basic assumptions of the model

The numerical simulation analysis of the model calculation process conforms to the following four assumptions [26]: (1) the soil in numerical models is saturated; the water in the soil seepage process conforms to the Darcy law; (2) soil grain and pore water is incompressible; (3) set the material property of the gravel piles as isotropic linear elasticity; (4) considering the numerical stability of the PL-Finn model, the compression coefficient and permeability coefficient of the soil are set as constant, and the permeability coefficients are equal in each direction.

3.5 The setting of monitoring points

To monitor structural displacement and the excess pore water pressure ratio in the medium-coarse sand layer, monitoring points are installed on both sides of the diversion dike top, while monitoring units for this ratio are embedded within the sand layer. Fig 7 is the Schematic diagram of monitoring points and units. Table 6 is model monitoring points and element coordinates.

Table 6: Model monitoring points and units coordinates.

Schematic diagram of monitoring points and units.

3.6 Analysis procedure

When simulating the dynamic response of a composite foundation under earthquake action, static analysis and seepage flow calculation of the model should be conducted first to obtain the correct initial ground stress field and seepage field. During dynamic calculation, issues such as material model and parameters, boundary conditions, and damping types should be taken into consideration. The specific analysis procedure flowchart is shown in Fig 8.

The analysis procedure.

4 Analysis of numerical simulation results

4.1 Analysis of the residual deformation

The residual deformation cloud charts for the conditions before and after gravel pile reinforcement are shown in Figs 9 and 10, with deformation time-history curves at monitoring points shown in Fig 11. Comparative analysis of residual deformation under the two working conditions indicates that liquefaction in bilateral sand layers induces symmetrical lateral displacement of the diversion dike about its central axis, as evidenced by the horizontal residual deformation cloud chart. In contrast, the underlying sand layer exhibits minimal displacement due to the absence of liquefaction, while maximum horizontal displacement occurs at both sides of the structure.

Cloud chart of residual deformation distribution before reinforcement (a) Cloud chart of horizontal residual deformation distribution (b) Cloud chart of Vertical residual deformation distribution.

Cloud chart of residual deformation distribution after reinforcement (a) Cloud chart of horizontal residual deformation distribution (b) Cloud chart of Vertical residual deformation distribution.

Time-history curve of deformation at monitoring points before and after reinforcement (a) Horizontal displacement deformation curve of point a (b) Vertical displacement deformation curve of point a (c) Horizontal displacement deformation curve of point b (d) Vertical displacement deformation curve of point b.

According to the time-history curve of horizontal displacement deformation of the monitoring points, the horizontal displacement of monitoring point a is 1.11m and that of monitoring point b is 1.0m before reinforcement. After reinforcement, the horizontal displacement of monitoring point a is 0.72m, and monitoring point b is 0.69m. The residual deformation of horizontal displacement is reduced by 35%.

From the vertical residual deformation distribution cloud chart, the whole diversion dike has uniform settlement under earthquakes, and the largest settlement occurs at the top of the diversion dike. The foundation settlement range below the structure gradually decreases with the increase of depth, and the heave deformation occurs on the contact position between the two sides of the structure and the mucky soil. It can be found that from the time-history curve of vertical displacement deformation, before reinforcement, the subsidence of the monitoring points a and b is 1.18m and 1.17m respectively. After reinforcement, the subsidence is reduced to 0.62m and 0.60m. The residual deformation of vertical settlement is reduced by 47%. According to the reduction of vertical and horizontal residual deformation of the diversion dike, it can be proposed that the seismic performance of the diversion dike has been significantly improved.

4.2 Analysis of the excess pore water pressure

When saturated sand is subjected to external loads, its volume decreases as the particles are compressed. This compression causes the pore water pressure to increase steadily, eventually exceeding the original hydrostatic pressure. The pressure in excess of the hydrostatic pressure is termed the excess pore water pressure. The development of excess pore water pressure is a key indicator of the liquefaction state in saturated sand.

Fig 12 shows the time-history curves of excess pore water pressure for each monitoring unit before and after reinforcement. Following the application of seismic waves at the model base, the excess pore water pressure at each unit rapidly peaks, then decreases slightly, and finally stabilizes. As the distance from the pile center increases, the excess pore water pressure rises, indicating a progressive weakening of the gravel piles’ drainage effect. Comparison shows that peak excess pore water pressures in the reinforced area are significantly lower after reinforcement, whereas peaks in the unreinforced area show little change. These results demonstrate that compacted gravel piles provide effective drainage, facilitating the dissipation of excess pore water pressure in the medium-coarse sand layer and thereby reducing liquefaction potential.

Time-history curves of excess pore water pressure for different monitoring units: (a) Element A (b) Element B (c) Element C (d) Element D (e) Element E.

4.3 Analysis of the excess pore water pressure ratio

In this study, FLAC^3D^ performs numerical simulations and analyses.The excess pore water pressure ratio, defined as the ratio of pore water pressure increment to lateral effective consolidation stress, serves as the liquefaction criterion to ensure computational accuracy. Experimental data [27] indicate that liquefaction initiation occurs when this ratio reaches 0.87 in silty soils, with complete liquefaction at 1.0. Herein, a threshold ratio of 0.8 is adopted to define the liquefaction state, implemented through FISH programming for concurrent computation.

Fig 13 compares the time-history curves of the excess pore water pressure ratio for each monitoring point before and after reinforcement. The results show that installing compacted gravel piles in the saturated sand foundation significantly reduces the peak excess pore water pressure ratio at monitoring points within the reinforced area, but causes little change outside it. Before reinforcement, the excess pore water pressure ratio rises rapidly at 7s and then falls rapidly, stabilizing at 14s. This indicates that sand liquefaction occurs (when effective stress reaches zero), followed by gradual dissipation of pore water pressure. During seismic loading, the excess pore water pressure ratio increases with increasing horizontal distance from the pile center, indicating a progressive weakening of the gravel piles’ drainage effect. Consequently, the anti-liquefaction performance of the gravel pile-reinforced composite foundation is significantly enhanced. Within the reinforced area, the peak excess pore water pressure ratio near the center is about 0.4, while near the edge it ranges between 0.45 and 0.5. Outside the reinforced area, the sand layer experiences varying degrees of liquefaction. Therefore, the goal of eliminating site liquefaction is largely achieved under the design-level earthquake intensity. Fig 14 shows the distribution of liquefied zones in the three-dimensional model for both (pre- and post-reinforcement) conditions.

Time-history curves of excess pore pressure ratio for different monitoring units: (a) Element A; (b) Element B; (c) Element C; (d) Element D; (e) Element E.

Distribution of liquefied zones in the three-dimensional model: (a) Before reinforcement; (b) After reinforcement.

5 Discussion

In this study, we adopt the PL-Finn constitutive model—a secondary development of the Finn model incorporating post-liquefaction flow deformation behavior. Numerical simulations of the diversion dike closely align with field investigation data. Compared to the original Finn model, the PL-Finn model more accurately captures deformation patterns and liquefaction zones in saturated sand under dynamic loading. All conventional mechanical parameters required for static and dynamic analyses can be determined through laboratory tests or empirical formulas, granting the model strong practical applicability in engineering.

Observations (Figs 9–10) reveal that post-liquefaction sand flow induces regular residual deformation: the structure undergoes near-symmetrical lateral sliding on both flanks and uniform upper-section settlement, with magnitude decreasing progressively with depth. Gravel pile reinforcement triggers lateral heave deformation in liquefiable soils outside the reinforced zone due to foundation subsidence-induced compression, exhibiting asymmetric patterns. These characteristics align with Zou’s findings [28]. Mechanistically, adjacent unreinforced soils exhibit limited drainage capacity, causing delayed dissipation of excess pore pressure post-liquefaction. Consequently, low volumetric strain and weak deformation resistance lead to significant lateral and upward displacement when compressed by overlying foundations.

Notably, Adampiro and Derakhshandi [29] investigated layered liquefiable soils’ effects on seismic site response via experiments and numerical simulations. Their findings clarify that the relative thickness of the liquefiable layer and overlying non-liquefied crust govern seismic wave attenuation, amplification, and ground motion intensity, consequently altering liquefaction characteristics. Fig 12 reveals higher peak excess pore pressure at the center versus sides, attributed to greater sand layer thickness—validating the thickness-liquefaction correlation. Liquefaction initiates when the excess pore pressure ratio exceeds 0.8 in sand layers under seismic loading (Fig 13). Given the unique requirements of nuclear power plant marine structures and the diversion dike’s critical role in operational safety, foundation reinforcement becomes paramount.

Additionally, Iranian scholars [30] demonstrated that gravel piles enhance drainage within about 2.5m from their center, with effects diminishing beyond this radius—corroborating our findings on effective liquefaction mitigation ranges. Their study recommends optimal settlement control at pile spacings of 2.5–3.5 times pile diameter [30], aligning with our design (Fig 4). Notably, El Sherbiny and Wang [31] determined square patterns outperform triangular or hexagonal arrangements for tunnel liquefaction mitigation. Closer pile spacing enhances soil stiffness and enables uniform drainage, thereby reducing excess pore pressure and tunnel uplift. While distinct from our focus, future work should assess pile arrangements (spacing, diameter, pattern) to optimize cost-effective liquefaction countermeasures.

Owing to the numerical stability constraints of the PL-Finn model, this study adopted the assumption of constant soil permeability. However, pore structure rearrangement during liquefaction may alter permeability, thereby affecting the dissipation of excess pore water pressure and the drainage efficiency of gravel piles. To address this knowledge gap, we will quantify the effects of permeability variations on pore water pressure dissipation and drainage efficiency in ongoing research (S1 File).

In this study, due to the conservative design philosophy inherent in nuclear power engineering, the pile-soil contact was modeled under ultimate bearing conditions assuming a perfectly smooth state with zero friction angle. This setting may lead to an underestimation of the lateral load transfer efficiency of stone columns, while potentially overestimating foundation settlement to some extent, resulting in conservative simulation outcomes. This approach aligns with the principle of prioritizing computational efficiency. Although the results obtained from the numerical simulations demonstrate significant potential for application in the design of nuclear marine infrastructure, these conclusions have not yet been validated through field data or benchmark testing due to current research limitations. To enhance the reliability and generalizability of the findings, we will conduct scaled shaking table model tests on diversion dikes in future research. Systematic comparison between measured dynamic responses and three-dimensional numerical simulation results will validate the rationality of constitutive model parameter selection and the effectiveness of the computational methodology.

6 Conclusion

Based on the actual liquefaction engineering of a coastal nuclear power project, this study studies the dynamic characteristics and anti-liquefaction reinforcement effect of the composite foundation of gravel piles in liquefiable sites, and carries out numerical analyses of the seismic dynamic response of the composite foundation-foundation-superstructure system of gravel piles in liquefiable sites. The main conclusions are as follows:

(1)Without reinforcement, liquefaction occurs; gravel pile reinforcement significantly enhances liquefaction resistance. In the reinforced zone, the maximum excess pore pressure ratio at the foundation base is about 0.50, preventing liquefaction under design-level ground motion.(2)During the earthquake, the excess pore water pressure at each monitoring point rapidly increases to the peak value, then decreases slightly, and finally hovers around the stable value. Gravel pile reinforcement significantly reduces the peak excess pore water pressure in the composite foundation, demonstrating that gravel piles can effectively dissipate excess pore water pressure and promptly discharge excess pore water.(3)While improving the liquefaction resistance of the reinforced area, gravel pile reinforcement greatly reduces the seismic settlement and displacement of liquefiable sites. Before reinforcement, the maximum lateral slip is about 1.11 m, and the maximum vertical settlement is about 1.18m. After reinforcement, the maximum lateral slip is about 0.72m, and the maximum vertical settlement is about 0.62 m. The vertical and horizontal residual deformation of the diversion dike is significantly reduced, and the stability of the structure is greatly improved.

Supporting information

S1 FilePL-Finn Model-based computational code for liquefaction analysis.(ZIP)

S1 FigPL-Finn model development procedure diagram.(TIF)