## 3 Numerical Model
### Finite Element Model  
The model was established using the finite element analysis software ABAQUS, as shown in Figure 11[32]. The calculation and analysis method adopted the overall dynamic time history analysis method of the soil-underground structure, based on equivalent linearization, as proposed by [PERSON] et al. [33]. The elastic modulus of the soil layer was calculated using Equations (1) and (2). \[G=\rho V_{S}^{2} \tag{1}\]  
\[E=2(1+\
u)G \tag{2}\] \begin{table}
\begin{tabular}{c c c c} \hline \hline
**Ground Level** & **Densities/kg/m\({}^{3}\)** & **Shear Wave Velocity/m/s** & **Poisson’s Ratio** \\ \hline sandy soil layer & 2000 & 380 & 0.3 \\ clay layer 1 & 2000 & 150 & 0.3 \\ clay layer 2 & 2000 & 125 & 0.3 \\ clay layer 3 & 2000 & 100 & 0.3 \\ clay layer 4 & 2000 & 75 & 0.3 \\ \hline \hline \end{tabular}
\end{table} 
Table 2: Physical properties of soils. Figure 10: Schematic diagram of the soil-underground structure. \begin{table}
\begin{tabular}{c c c c} \hline \hline
**Working Condition** & **Interlayer Thickness/m** & \
\begin{tabular}{c} **Shear Wave Velocity in Weak** \\ **Interlayers/m/s** \\ \end{tabular} &
\begin{tabular}{c} **Homogeneous Layer Shear** \\ **Wave Velocity/m/s** \\ \end{tabular} \\ \hline
0 & interlayer-free & / & 380 \\
1 & 1 & 125 & 380 \\
2 & 2 & 125 & 380 \\
3 & 3 & 125 & 380 \\
4 & 4 & 125 & 380 \\
5 & 3 & 150 & 380 \\
6 & 3 & 100 & 380 \\
7 & 3 & 75 & 380 \\ \hline \hline \end{tabular}
\end{table} 
Table 3: List of operating conditions. Here, \(\rho\) represents the density of the soil layer, \(V_{\rm s}\) stands for the shear wave velocity of the soil layer, and \(\
u\) represents the Poisson's ratio of the soil layer. In defining the material properties of the soil layer, the seismic mitigation of the soil material needed to be specified. In this model, the Rayleigh seismic mitigation model was selected. The frequency analysis module in ABAQUS was used to obtain the natural frequency of the test soil. The first-order and second-order frequencies were selected, and the Rayleigh seismic mitigation coefficient was calculated using Formulas (3) and (4). \[\alpha=\frac{2D\times w_{1}\times w_{2}}{w_{1}+w_{2}} \tag{3}\]  
\[\beta=\frac{2D}{w_{1}+w_{2}} \tag{4}\]  
Here, \(w_{1}\) and \(w_{2}\) are the first and second order self-oscillation frequencies of the soil layer and \(D\) is the seismic mitigation ratio taken as 0.1 [34]. The model materials were all elastic constitutive models. The total width of the soil-structure system model was 140 m, and the height was 40 m, which met the requirement of being more than seven times the width of the structure. The middle columns of the station were distributed equidistantly along the longitudinal direction of the station, with a spacing of 6 m between adjacent middle columns. Due to the two-dimensional model adopted in this paper, the station structure used the plane strain unit to establish the model, and the material of the middle column needed to be reduced according to Formula (5) [35]. \[E^{\prime}=E\frac{h}{h^{\prime}}\ \ \rho^{\prime}=\rho\frac{h}{h^{\prime}} \tag{5}\]  
Here, \(E^{\prime}\) is the equivalent elastic modulus, \(\rho^{\prime}\) is the equivalent material density, \(h\) is the length of the middle column section along the longitudinal direction of the subway, and \(h^{\prime}\) is the distance between the two adjacent middle columns in the longitudinal direction. After reduction, the elastic modulus of the middle column material is 4.67 GPa, the density is 333.3 kg/m\({}^{3}\), and the Poisson's ratio is 0.2. The binding boundary utilizes the MPC (Multi-Point Constraint) node degree of freedom coupling feature within ABAQUS, effectively binding the finite element model nodes of the soil to the boundary nodes at the same height, ensuring these nodes share identical displacements throughout the analysis [36]. Previous studies have shown that the boundary form has good calculation accuracy in soil-structure dynamic analysis [37] and pseudo static analysis [38].