Multi-Factor Inﬂuence Analysis on the Liquefaction Mitigation of Stone Columns Composite Foundation

: To optimize the design of stone columns composite foundation for liqueﬁable ground improvement in the Tibar Bay Port Project, a 3D Finite Element (FE) analysis is implemented on the earthquake response and liquefaction mitigation effect. Nine improvement schemes are designed with the orthogonal design method. Taking peak ground acceleration and peak excess pore pressure ratio as the target indicators, the inﬂuences of four factors, including diameter, replacement ratio, stiffness, permeability ratio, of stone columns are analyzed by means of range analysis, and subsequently, the optimal ground improvement design is obtained. The analysis results indicate that the responses of ground acceleration and excess pore pressure ratio are relatively sensitive to stone columns’ permeability ratio and a little sensitive to the replacement ratio. The stiffness and diameter ranging in the prescribed boundary only have negligible effect. The mitigation effect of drainage is rather signiﬁcant when the ratio of the stone columns’ permeability to the soils’ permeability is greater than 100.


Introduction
Ground liquefaction has a destructive effect on infrastructure, and disasters caused by soil liquefaction during earthquakes take place frequently around the world [1][2][3][4]. At present, the commonly used ground improvement measures in engineering practice for liquefaction mitigation include densification, refilling, and granular columns among others [5]. Since Seed and Booker [6] first studied the drainage effect of granular columns on liquefaction mitigation in 1977, the ground improvement method of granular columns has been gradually adopted in engineering projects around the world [7,8]. Extensive proofs have been obtained from previous research including the earthquake case histories and experimental investigations for the effectiveness of granular columns on liquefaction mitigation [8][9][10]. Nowadays it is recognized that the mitigation mechanism of granular columns consists of three aspects: (1) densification, (2) drainage, and (3) shear reinforcement [10][11][12][13]. The installation of granular columns will densify surrounding soils. Furthermore, their relatively higher modulus and stiffness will be able to help to confine the ground deformation as well as reduce the proportion of earthquake-induced shear stress in the surrounding soils. Moreover, the high permeability of granular columns further facilitates the dissipation of excess porewater pressure (EXPP). Despite the previous efforts in theoretical and experimental research, the standard design approach for granular column composite foundation has not been proposed. As made up of granular material, the granular columns' mechanical behaviors and deformation patterns are closely related to the mechanical properties and deformation of surrounding soils, which adds complexity to investigating the earthquake response characteristics of the composite foundation and generalizing an applicable design approach. Therefore, pertinent studies are still needed for granular columns in particular engineering sites.
The dynamic numerical analysis is a fundamental means to study the earthquake response of granular columns composite foundation in liquefiable soils, among which the FE analysis based on Biot's dynamic consolidation equation can fully couple the pore pressure with the dynamic response of the soil skeleton and becomes one of the most extensively used approaches [14]. In order to assess the effectiveness and mitigation of stone columns, three-dimensional FE simulations were implemented on OpenSees in various studies [15][16][17][18][19] with the same computational formulation. These studies assumed the saturated soils to be solid-fluid material following the u-p formulation which is developed based on Biot's theory by Chan [20] and Zienkiewicz et al. [21]. The u denotes displacement of the soil skeleton and p denotes pore pressure. The multiple yield surface elastoplastic constitutive model, PDMY (PressureDependMultiYield), was adopted in these studies to simulate the materials of both surrounding soils and granular columns. The seismic responses of the granular column composite foundation were simulated. Furthermore, the liquefaction mitigation effects were analyzed such as the shear reinforcement, and the lateral displacements restraining. Tang et al. [22] applied the same numerical method to simulate the centrifuge test of the saturated silty site improved with granular columns and analyzed the dynamic characteristics such as ground acceleration and pore pressure. Hereby, the results indicated that the numerical simulation was in good agreement with the experimental measurement. There are also FE simulations based on u-p formulations with other constitutive models, such as the hypoplastic model [23] and critical state model based on generalized plasticity theory [24], which are conducted to evaluate the performance of granular columns in liquefaction mitigation. Additionally, Zou et al. [25] used FLAC3D for dynamic response analysis, applied the elastoplastic constitutive model which is able to predict the large post liquefaction deformation and equivalent nonlinear model to simulate materials of sand and granular columns, respectively, and evaluated the liquefaction mitigation effect in terms of shear reinforcement and drainage of the granular columns.
The current study is based on the background of design optimization for stone columns composite foundation at a liquefiable site in the East Timor Tibar Bay port project (TBPP). The project, 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 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, the gradations of which are spatially varied, and a high proportion of fine sands can be found in local areas. To mitigate the potential liquefaction, granular columns are adopted for the ground improvement and are mainly installed in areas where layers of sand or sand-gravel mixture a very loose to loose or partially dense state spread. The standard penetration test (SPT) values of these layers range from 0 to 18, with a mean value being 7. The design relative density of the site is 50%, and the design permeability is 6.4 × 10 −2 cm/s.
To investigate the influence of relevant factors on mitigation and further optimize the design, the stone columns are firstly designed using the orthogonal design method. Then, seismic response analysis is conducted through fully coupled FE simulations on OpenSees, which adopt the aforementioned u-p formulation for saturated soils. The constitutive model of soils plays a critical role in this fully coupled numerical approach and decides its prediction capability. The multi-yield surface elastoplastic constitutive 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 [15][16][17][18][19]26,27]. In addition, with a total of 17 input parameters, the PDMY model allows for a comprehensive multi-factor analysis. Therefore, PDMY was adopted to simulate the liquefiable sand and stone materials. Taking seismic response of ground acceleration and excess pore pressure ratio (EXPPR) as key indexes, the influences of design parameters including columns diameter, replacement ratio, shear stiffness, and permeability ratio on mitigation effect are evaluated. Subsequently, the design scheme of the stone columns is optimized based on the analysis results. Overall, the conducted research provides a means of comparing the significance of relevant factors to the mitigation effect of stone columns, and the insights and conclusions can be helpful for practical engineering design.

Design Schemes of Stone Columns
In order to obtain the optimal scheme of the stone columns, the orthogonal method was applied for the scheme design, and numerical analyses of multi-factor influences were carried out. For discrete single objective optimization problem with m factors varying in k levels, a total of m k combinations of levels would exist in the search space. To address the problem when m and k are large and analyzing every single combination is rather costly, the orthogonal experimental design method is developed. Through an orthogonal table, the method can provide selected combinations which are uniformly scattered over the space of all possible combinations [28]. The orthogonal table is a fractional factorial array, where each row represents each combination of factors in corresponding levels, and each column represents the relevant factor that is investigated and can be changed within the range of levels in the experiments. One column can be evaluated independently from the others, for which it is named as an orthogonal table. The results obtained from each row are further analyzed to evaluate the significance of factors and estimate the optimal combination of levels of each factor.
Four factors of the stone columns are selected as the major variables in the numerical experiments: column diameter, replacement ratio, shear stiffness, and column permeability ratio. Each factor varied on three levels, the boundary values of which are determined by actual engineering practice, see Table 1. The L9 (3 4 ) orthogonal table is designed, in which four columns are designated for investigated factors, and each factor ranges corresponding to the prescribed three levels. Based on that, nine schemes of different combinations of factors are obtained. A set of free field models was considered to provide a benchmark for comparing the liquefaction mitigation effect among the groups of stone columns. The numerical experiment program was finally formed as listed in Table 2.

Constitutive Model
PDMY02, which is developed based on the multisurface-plasticity theory for frictional cohesionless soils proposed by Prevost [29][30][31], is applied for modeling both sand and stone materials. It defines a conical yield surface as presented. The outermost yield surface is the failure surface, and several similar yield surfaces exist inside. During the loading, the stress state moves according to the corresponding hardening law and flow rules. A combination of associative and non-associative flow rules is adopted for PDMY02, which decomposes the plastic potential surface and the outer normal to the yield surface into deviatoric components and volumetric components. The deviatoric components follow the associative flow rules whereas the volumetric components follow the non-associative flow rules. Non-associativity is proposed to describe the coupling relationship between the shear stress and the volume strain during cyclic loading. In order to generate the hysteretic response of soil under cyclic shear loads, the model obeys the kinematic hardening law proposed by Mroz [32] and makes remediation to enhance the computational efficiency as well as model robustness [33]. Table 2. Scheme of simulation models.

No. Diameter Replacement Ratio Shear Stiffness Permeability Ratio
* The numbers in the table denote the different levels for each factor.
The PDMY02 model has a total of 17 input parameters. Yang et al. [34] calibrated the parameters of sand with different densities, which are referred to in this paper for the sand materials. The measured data at the in situ site are also utilized to obtain the parameters. The parameters of the material of the stone columns are determined using the calibration method adopted by Rayamajhi et al. [17] as well as by the design scheme proposed in the last section, as listed in Table 3.  Noted that the elastic shear modulus of the stone columns is determined by the design schemes proposed in the last section; the elastic bulk modulus can be calculated as follows [35], where G r is the elastic shear modulus, and ϑ is the Poisson's ratio which is taken as 0.33.

Finite Element Model
The 3D dynamic analysis was carried out using OpenSees. The stone columns are arranged in a square array, which is highly periodic and symmetric in structure, so the numerical analysis can be conducted with the model of one-half unit cell [15][16][17][18], as illustrated in Figure 1. The model is overall 30 m high. The top layer of 5 m thickness is the dense saturated sand with 80% relative density, and the underlain layer of 25 m is medium dense saturated sand with 50% relative density. This layout is aimed to simulate the in situ Appl. Sci. 2022, 12, 7308 5 of 11 soil profile with the dense reclamation layer on top of the natural liquefiable sand. The groundwater table is set at the surface level. The permeability k sand of all sand materials is 6.4 × 10 −2 cm/s. The length of the columns is 15 m, and the replacement ratio is defined as the ratio of the cross-sectional area of the column to the unit cell. The center-to-center distance of columns for each scheme is calculated according to the replacement ratio and the column diameter. Other parameters of sand and stone materials are documented in Tables 1 and 3. The site is modeled using brickUP elements [34], following the u-p dynamic consolidation equation, while the thickness of elements is 1 m.

Finite Element Model
The 3D dynamic analysis was carried out using OpenSees. The stone columns are arranged in a square array, which is highly periodic and symmetric in structure, so the numerical analysis can be conducted with the model of one-half unit cell [15][16][17][18], as illustrated in Figure 1. The model is overall 30 m high. The top layer of 5 m thickness is the dense saturated sand with 80% relative density, and the underlain layer of 25 m is medium dense saturated sand with 50% relative density. This layout is aimed to simulate the in situ soil profile with the dense reclamation layer on top of the natural liquefiable sand. The groundwater table is set at the surface level. The permeability of all sand materials is 6.4 × 10 −2 cm/s. The length of the columns is 15 m, and the replacement ratio is defined as the ratio of the cross-sectional area of the column to the unit cell. The center-to-center distance of columns for each scheme is calculated according to the replacement ratio and the column diameter. Other parameters of sand and stone materials are documented in Tables 1 and 3. The site is modeled using brickUP elements [34], following the -dynamic consolidation equation, while the thickness of elements is 1 m. Considering the periodicity of the model [36], 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 displacement of nodes on the lateral boundaries are fixed in -direction; (3) the pore pressure on the soil surface is fixed to be zero; and (4) the nodes at the base are fixed, and the input seismic loading is imposed on the base along the -direction. The earthquake history recorded by the near-site seismic station is selected as the input seismic loading, and the peak ground acceleration (PGA) is adjusted to 0.6 g in order to meet the requirements corresponding to the seismic fortification intensive. The acceleration time history curve is shown in Figure 2. Considering the periodicity of the model [36], 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 z-direction; (2) the displacement of nodes on the lateral boundaries are fixed in y-direction; (3) the pore pressure on the soil surface is fixed to be zero; and (4) the nodes at the base are fixed, and the input seismic loading is imposed on the base along the x-direction. The earthquake history recorded by the near-site seismic station is selected as the input seismic loading, and the peak ground acceleration (PGA) is adjusted to 0.6 g in order to meet the requirements corresponding to the seismic fortification intensive. The acceleration time history curve is shown in Figure 2.

Computation Procedure
Apply uniform linear elastic mechanical properties 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, transform the corresponding input parameters of sand and stone materials and perform elastoplastic analysis to obtain the initial stress field and pore pressure field. The ground motion is imposed on the base of the model by acceleration

Computation Procedure
Apply uniform linear elastic mechanical properties to all the materials in the model to establish the initial static stress condition as well as to avoid the stress concentration Appl. Sci. 2022, 12, 7308 6 of 11 or arching effects. After that, transform the corresponding input parameters of sand and stone materials and perform elastoplastic analysis to obtain the initial stress field and pore pressure field. The ground motion is imposed on the base of the model by 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 conditions.

Simulation Results
The earthquake responses of different schemes of stone columns are analyzed. The ground acceleration of point 1, indicated in Figure 1c

Computation Procedure
Apply uniform linear elastic mechanical properties to all the materials in the to establish the initial static stress condition as well as to avoid the stress concentra arching effects. After that, transform the corresponding input parameters of san stone materials and perform elastoplastic analysis to obtain the initial stress field an pressure field. The ground motion is imposed on the base of the model by accel time history, namely the vibration method. During the computation, Rayleigh da of 0.5% was applied to control high-frequency noise and ensure the damping ener sipation of the soils under small strain conditions.

Simulation Results
The earthquake responses of different schemes of stone columns are analyze ground acceleration of point 1, indicated in Figure 1c   For the free field model (FF-0), upon seismic loading, the EXPPR at depths of 2 m, 7 m, and 12 m climb up immediately and the peak values were reached in around 10 s, after which the pore water pressure was difficult to dissipate. The peak value of the EXPPR at each depth is approximate to 1.0, indicating that the sand has nearly fully liquefied. However, the peak value of the EXPPR at the depth of 2 m is slightly lower than that of deeper layers, which can be explained by the dense sand layer's stronger liquefaction resistance.
A reduction in peak values of the EXPPR at the depths of 2 m, 7 m, and 12 m can be seen in all three stone column models (SC-1/2/3). In addition, the presence of stone columns accelerates the dissipation of the pore pressure that starts to drop right after the peak amplitude of the seismic excitation. There is a relatively minor difference in the EXPPRs at the three different depths for the same model. However, the reduction in EXPPR of the SC-1 is apparently smaller than that of the SC-2 and 3. In terms of ground acceleration response, for the free field model, the amplitude of the ground acceleration decreases rapidly in 10 s, representing a characteristic of liquefaction, the vibration isolation. With regard to the model SC-1/2/3, the "vibration isolation" effect is somewhat suppressed, but the PGAs go up to a certain degree. The SC-2 and 3 exhibit a similar trend with the free field model during the beginning of 8 s concerning the ground acceleration. A similar response of SC-1 with free field lasts slightly longer, for about 10 s. The earthquake response of models SC-4 to SC-9 is not analyzed in detail hereby. The PGAs at the same observation points of each model and also the peak value of the EXPPR at different depths are collected in Table 4. For the free field model (FF-0), upon seismic loading, the EXPPR at depths of 2 m, 7 m, and 12 m climb up immediately and the peak values were reached in around 10 s, after which the pore water pressure was difficult to dissipate. The peak value of the EXPPR at each depth is approximate to 1.0, indicating that the sand has nearly fully liquefied. However, the peak value of the EXPPR at the depth of 2 m is slightly lower than that of deeper layers, which can be explained by the dense sand layer's stronger liquefaction resistance.
A reduction in peak values of the EXPPR at the depths of 2 m, 7 m, and 12 m can be seen in all three stone column models (SC-1/2/3). In addition, the presence of stone columns accelerates the dissipation of the pore pressure that starts to drop right after the peak amplitude of the seismic excitation. There is a relatively minor difference in the EXPPRs at the three different depths for the same model. However, the reduction in EXPPR of the SC-1 is apparently smaller than that of the SC-2 and 3. In terms of ground acceleration response, for the free field model, the amplitude of the ground acceleration decreases rapidly in 10 s, representing a characteristic of liquefaction, the vibration isolation. With regard to the model SC-1/2/3, the "vibration isolation" effect is somewhat suppressed, but the PGAs go up to a certain degree. The SC-2 and 3 exhibit a similar trend with the free field model during the beginning of 8 s concerning the ground acceleration. A similar response of SC-1 with free field lasts slightly longer, for about 10 s. The earthquake response of models SC-4 to SC-9 is not analyzed in detail hereby. The PGAs at the same observation points of each model and also the peak value of the EXPPR at different depths are collected in Table 4.

Influence Analysis and Scheme Optimization
Based on the computation results in Table 4, the sensitivity analysis was carried out by the range analysis method. The influences of the stone columns' design parameters on the peak values of the ground acceleration and the EXPPR in the surrounding soils were evaluated. Three mean values corresponding to three levels were calculated firstly for each factor on the basis of a particular set of computation results. The set is defined to be the group of results calculated from the schemes that have the same level for a certain factor, e.g., the three mean values for diameter corresponding to levels 1-3 are calculated from the computation results of SC-1/2/3, SC-4/5/6, SC-7/8/9, respectively. Then, the maximum difference (MD) of mean values for each factor is obtained, which could reflect the impact of the factor on the target index. Table 5 is the range analysis table of ground acceleration. It can be concluded that the influence of the permeability is the most significant, followed by the stone column's replacement ratio, while the diameter and shear stiffness of the stone column have little effect on the ground acceleration. When the permeability ratio of the stone columns rises from 1 cm/s (L1) to 5 cm/s (L2) and 10 cm/s (L3), the corresponding PGA goes up significantly, from 1.96 m/s 2 to 2.64 m/s 2 and 2.86 m/s 2 , respectively. When the replacement ratio increased from 10.1% (L1) to 12.9% (L2) and 15.7% (L3), the PGA climbed up slightly from 2.40 m/s 2 to 2.48 m/s 2 and 2.58 m/s 2 , respectively. The sensitivity of the EXPPR to various factors at different depths of 2 m, 7 m, and 12 m presents the same characteristics, for which only the range analysis results at the depth of 7 m are displayed hereby (see Table 6). The table indicates that the EXPPR is the most sensitive to the permeability ratio, which is followed by the replacement ratio, and the increase of the column diameter and shear stiffness can hardly reduce the peak value of EXPPR. When the permeability ratio of the stone columns rises from 1 cm/s (L1) to 5 cm/s (L2) and 10 cm/s (L3), the corresponding mean values of peaks of EXPPR drop from 0.85 to 0.35 and 0.22, respectively. A conspicuous decline in the EXPPR can be seen when the permeability ratio increases from 1 cm/s (L1) to 5 cm/s (L2), whereas a relatively small decline appears when permeability ratio continues increasing from of 5 cm/s (L2) to 10 cm/s (L3). With regard to the impact of the replacement ratio, when it increases from 10.1% (L1) to 12.9% (L2) and 15.7% (L3), the mean values of peaks of EXPPR decrease slightly from 0.53 to 0.47 and 0.43, respectively. The design of stone columns is optimized to achieve the objective of enhancing the mitigation effect while reducing the economic cost. Therefore, the EXPPR, which indicates the risk of liquefaction, is firstly taken as the major index. From the range analysis above, it can be proved that the permeability ratio of the stone columns is the most sensitive factor for the EXPPR of the site. When the permeability ratio of the stone columns ascended from 1 cm/s to 5 cm/s, the peak value of the EXPPR in the surrounding soils descended from 0.85 to 0.35, and the liquefaction risk has been greatly reduced. Furthermore, under the circumstance of a high permeability ratio, being 5 cm/s, the change of other design parameters has a relatively small effect on the reduction of the EXPPR, besides, continuing to increase the permeability ratio brings about a small decrease in the peak EXPPR. Therefore, the optimal scheme for the stone columns composite foundation of the project site is that the column diameter is 0.8 m, the replacement ratio is 10.1%, the shear stiffness is 2.20 × 10 5 kPa, and the permeability coefficient is 5 cm/s.

Recommendations for Practice
The conducted research analyzed the impact of factors on mitigation effect and made comparisons among the factors with regard to the sensitivity, which provided insights for improving the design of stone columns in practice. The analysis framework is also valuable and referable in relevant investigations.
The simulations present that stone columns were effective in reducing the peak values of EXPPR and accelerating the dissipation of pore pressure in surrounding soils. The stone column's permeability ratio was found to be a critical factor in this mitigation effect. The increase in replacement ratio can also lead to a reduction in the EXPPRs, but it is less significant than the permeability ratio. In contrast, the change of diameter and shear stiffness in the selected boundary exert a trivial impact on the peak values of EXPPRs. Similar conclusions were obtained in terms of the PGAs. The permeability ratio plays the most significant role in affecting the PGAs, which is followed by the replacement ratio, while the influence of diameter and shear stiffness is negligible. Additionally, the results regarding PGA manifested a pronounced difference against those regarding EXPPR, that the increase of permeability ratio and replacement ratio would result in larger PGAs. Such phenomena can be explained by the suppressing of vibration isolation due to the strengthened liquefaction resistance of surrounding soils.
The obtained results suggest that to mitigate the liquefaction risk by reducing the EXPPR, it is vital for stone columns to have a reasonable permeability ratio and replacement ratio. The increase of permeability is more effective to reduce the risk. Furthermore, the design can be further improved by adjusting the diameter as well as shear stiffness for economic reasons. However, it should be noted that the recommendations are made on the basis of the soils with particular mechanical properties in the investigated engineering case. Furthermore, because the influence of the factors is studied by changing the parameters in a limited range, more analyses are needed when the design parameters exceed the current investigated domain.
Overall, the conducted research does not consider the effects of construction processes on the surrounding soils (e.g., densification), for which the positive impact of replacement ratio may be underestimated and more efforts can be made in this aspect for further research. Besides, the drainage capacity of stone columns, which are suggested to be critical for liquefaction mitigation in the current study, could be undermined in the long run due to multiple effects such as clogging. As not addressed hereby, more attention should be given to these effects to properly evaluate the long-term serviceability of stone columns and improve the design and construction in practice.

Conclusions
In the current study, a 3D FE fluid-solid fully coupled dynamic analysis is conducted on seismic response of the stone columns composite foundations in liquefiable soils. The orthogonal method is applied to design the stone columns and work out the numerical experiment program firstly. Then, range analysis is adopted to evaluate the influence of factors on the mitigation effect of the stone columns. The design of the stone columns is optimized subsequently. The major conclusions are drawn as follows: (1) The stone columns can effectively reduce the peak value of the EXPPR in the surrounding soils and accelerate the pore pressure dissipation. In each simulation case, the time histories of EXPPR within the depth of stone columns present a similar development pattern and have approaching peak values.
(2) The presence of stone columns will however lead to a larger PGA compared with that in free field, which can be attributed to the suppressing of vibration isolation and therefore higher transmitted ground accelerations. As a result, it could not be always beneficial to install stone columns to improve liquefiable ground which may underlie upper structures. (3) The peak values of EXPPR and the PGA are most sensitive to the permeability ratio of the stone columns among the investigated factors, which is followed by the replacement ratio, while the column diameter and shear stiffness ranging in the prescribed boundary cause negligible influence. (4) The reduction in peak values of EXPPRs and the rise in PGAs caused by the permeability ratio are rather significant when the ratio of stone columns' permeability to surrounding soils' reaches about 100 times. The ongoing increase of the ratio after that brings about a relatively smaller positive effect.