## 2 Numerical Model
### 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  
\begin{table}
\begin{tabular}{c c c c c} \hline \hline
**No.** & **Diameter** & **Replacement Ratio** & **Shear Stiffness** & **Permeability Ratio** \\ \hline FF-0 & - & - & - & - \\ SC-1 & 1 \({}^{*}\) & 1 & 1 & 1 \\ SC-2 & 1 & 2 & 2 & 2 \\ SC-3 & 1 & 3 & 3 & 3 \\ SC-4 & 2 & 1 & 2 & 3 \\ SC-5 & 2 & 2 & 3 & 1 \\ SC-6 & 2 & 3 & 1 & 2 \\ SC-7 & 3 & 1 & 3 & 2 \\ SC-8 & 3 & 2 & 1 & 3 \\ SC-9 & 3 & 3 & 2 & 1 \\ \hline \hline \end{tabular} \({}^{*}\)The numbers in the table denote the different levels for each factor.  
\end{table}
Table 2: Scheme of simulation models.  
\begin{table}
\begin{tabular}{c c c c} \hline \hline \multirow{2}{*}{**Parameters**} & \multicolumn{2}{c}{**Medium Dense**} & **Dense Sand** & \multirow{2}{*}{**Stone**} \\  & **Sand** & **Dr = 80\%** & & \\ \hline Density/ (t\(\cdot\)m\({}^{-3}\)) & 1.9 & 2.1 & 2.14 \\ Void ratio & 0.7 & 0.55 & 0.45 \\ Shear modulus/kPa & \(1.00\times 10^{5}\) & \(1.30\times 10^{5}\) & - \\ Bulk modulus/kPa & \(2.33\times 10^{5}\) & \(2.60\times 10^{5}\) & - \\ Friction angle/(\({}^{\circ}\)) & 33.5 & 36.5 & 48 \\ PT angle/(\({}^{\circ}\)) & 25.5 & 26 & 30 \\ Peak shear strain & 0.1 & 0.1 & 0.1 \\ Reference confining pressure/kPa & 101 & 101 & 101 \\ Pressure dependent coefficient & 0.5 & 0.5 & 0.5 \\ Contraction parameter & 0.045 /5.0/0.15 & 0.013/5.0/0.0 & 0.005/0.5/0.0 \\ Dilation parameter & 0.06/3.0/0.15 & 0.33/0.0/0.0 & 0.40/3.0/0.0 \\ Liquefaction parameter & 1.0/0 & 1.0/0 & 1.0/0 \\ \hline \hline \end{tabular}
\end{table}
Table 3: Parameters of PDMY02 material.  
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 \(\times\) 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.  
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.