Influence Analysis of Liquefiable Interlayer on Seismic Response of Underground Station Structure

To study the influence law of the seismic response of underground station structures at liquifiable interlayer sites, a two-dimensional numerical model of the interaction between the soil and station structure was established based on the finite difference software FLAC3D. The nonlinear dynamic response of the station structure located at the liquifiable interlayer site was analyzed considering the location distribution, relative density, and thickness of the liquifiable interlayer. The results show that the deformation of the structure is greatest when the liquifiable interlayer is distributed on both sides of the station side walls, while the interlayer has an energy-dissipating and damping effect on the upper station structure when it is located at the bottom of the structure. The lower the relative density of the liquifiable interlayer is, the stronger the internal dynamic response of the structure will be, and the more unfavorable it will be to the seismic resistance of the structure. When the liquefiable interlayer is only present in the lateral foundation of the station, an increase in its thickness results in a stronger shear effect on the structure and a higher probability of damage. However, when the thickness of the liquifiable interlayer reaches a point where the entire station is placed within it, the lateral force and deformation of the structure are significantly reduced.


Introduction
Large-scale transportation construction will inevitably result in underground space structures crossing liquefiable soils. When an earthquake occurs, the liquefiable soil will "liquefy" under the vibration load, causing the site soil to lose its bearing capacity and assume a liquid state. The liquefaction of the site soil will have an important impact on the underground structure [1,2]. For example, the 1995 Hanshin earthquake in Japan caused the devastating collapse of the Daikai subway station. The investigation found that extreme liquefaction occurred at both the coastal backfill site and in the artificial island area, causing widespread, permanent deformation of the foundation in the horizontal direction [3][4][5]. In the 2008 Wenchuan earthquake in Sichuan, the phenomenon of foundation liquefaction occurred in specific areas with severe earthquake damage, and underground structures such as tunnels were damaged to varying degrees [6,7]. In the 2011 earthquake off the Pacific Ocean in Japan, soil liquefaction occurred in a large region of coastal areas, many buildings were subjected to serious foundation subsidence, and some shallowly buried underground structures even appeared to rise above the ground [8,9]. Therefore, it is of great practical importance to study the seismic performance of underground structures at liquefiable sites.
In recent years, scholars in related fields have successively studied the seismic response law of underground station structures in liquefiable strata. Chen et al. [10] conducted a series of shaking table experimental studies on irregular cross-section subway underground structures at liquefiable sites and provided a specific description of the liquefaction development law of foundations and the damage mechanism of station structures. Zhuang Appl. Sci. 2023, 13, 9210 2 of 17 et al. [11] revealed the seismic response characteristics of underground station structures in slightly inclined liquefiable foundations through shaking table tests and further compared and analyzed the effect of ground inclination changes on station structures in liquefiable sites. Ding et al. [12] conducted a series of shaking table tests to reveal the effect of groundwater level changes on the degree of soil liquefaction and the seismic response pattern of underground structures. In addition to the methods employed through model tests, many scholars have used numerical simulations to investigate the dynamic response characteristics of subsurface structures at liquefiable sites [13][14][15][16][17][18]. Bao et al. [13] interpreted the seismic performance and uplift mechanism of large subway station structures in liquefiable soils considering the excavation process during the construction of the structure. Hu et al. [14,15] used the numerical analysis method coupled with the finite element and finite difference method to analyze the dynamics of rectangular structures in saturated sandy soil sites and explored the influences of different relative densities of sandy soils on the seismic response and uplift effect of underground structures. Tian et al. [16] explored the force and deformation mechanism of an underground corridor structure at a non-homogeneous liquefiable site and further investigated the effect of burial depth on the dynamic response of the structure.
In the above studies, the foundations were mostly considered as homogeneous liquefiable sites, or the station structures were wholly placed in liquefiable layers. The seismic response characteristics of station structures in non-homogeneous layered liquefiable sites need to be investigated further. Shen et al. [17] and Chen et al. [18] showed that underground structures undergo much more seismic damage at laminated liquefiable sites than at homogeneous liquefiable sites. In engineering practice, liquefiable foundations are mostly present in the form of liquefiable interlayers. The existing research is not explicit with respect to the influence law of liquefiable interlayers on the seismic response of underground station structures. In view of this, we take the liquefiable interlayer site as our research object, integrate the influences of the spatial location distribution, relative density and thickness of the liquefiable interlayer on the structure of a subway station under earthquake action, and draw some practical conclusions for seismic stability analysis in similar projects.

Geometric Model and Site Condition Setting
A typical subway station is a two-story, three-span structure. The cross-sectional dimensions of the main part are shown in Figure 1. The top plate of the station structure is buried at a depth of 5 m, the width of the structural model is 21.34 m, the height is 14.89 m, the size of the central column of the station structure is 0.6 m × 1.1 m, and the distance between the longitudinal columns of the central column is 9.12 m. uefaction development law of foundations and the damage mechanism of station struc tures. Zhuang et al. [11] revealed the seismic response characteristics of underground sta tion structures in slightly inclined liquefiable foundations through shaking table tests and further compared and analyzed the effect of ground inclination changes on station struc tures in liquefiable sites. Ding et al. [12] conducted a series of shaking table tests to revea the effect of groundwater level changes on the degree of soil liquefaction and the seismic response pattern of underground structures. In addition to the methods employed through model tests, many scholars have used numerical simulations to investigate the dynamic response characteristics of subsurface structures at liquefiable sites [13][14][15][16][17][18]. Bao et al. [13] interpreted the seismic performance and uplift mechanism of large subway sta tion structures in liquefiable soils considering the excavation process during the construc tion of the structure. Hu et al. [14,15] used the numerical analysis method coupled with the finite element and finite difference method to analyze the dynamics of rectangular structures in saturated sandy soil sites and explored the influences of different relative densities of sandy soils on the seismic response and uplift effect of underground struc tures. Tian et al. [16] explored the force and deformation mechanism of an underground corridor structure at a non-homogeneous liquefiable site and further investigated the ef fect of burial depth on the dynamic response of the structure.
In the above studies, the foundations were mostly considered as homogeneous liq uefiable sites, or the station structures were wholly placed in liquefiable layers. The seis mic response characteristics of station structures in non-homogeneous layered liquefiable sites need to be investigated further. Shen et al. [17] and Chen et al. [18] showed that un derground structures undergo much more seismic damage at laminated liquefiable sites than at homogeneous liquefiable sites. In engineering practice, liquefiable foundations are mostly present in the form of liquefiable interlayers. The existing research is not explici with respect to the influence law of liquefiable interlayers on the seismic response of un derground station structures. In view of this, we take the liquefiable interlayer site as our research object, integrate the influences of the spatial location distribution, relative density and thickness of the liquefiable interlayer on the structure of a subway station under earth quake action, and draw some practical conclusions for seismic stability analysis in similar projects.

Geometric Model and Site Condition Setting
A typical subway station is a two-story, three-span structure. The cross-sectional di mensions of the main part are shown in Figure 1. The top plate of the station structure is buried at a depth of 5 m, the width of the structural model is 21.34 m, the height is 14.89 m, the size of the central column of the station structure is 0.6 m × 1.1 m, and the distance between the longitudinal columns of the central column is 9.12 m. In this study, seven different site conditions were constructed, and three sets of numerical models were set up for comparison and analysis according to different variables. The influence of the location of the liquefiable interlayer on the dynamic response of the structure was investigated in cases a, b, and c; the influence of the relative density of the interlayer was investigated in cases b, d, and e; and the influence of the thickness of the interlayer was investigated in cases b, f, and g. The water level is assumed to be at the ground surface in each model, and the liquefiable interlayer is a sandy soil, while the rest of the site soil is non-liquefiable clay. A simplified model of the soil-station structure and the arrangement of the monitoring points are shown in Figure 2. The variable D is the burial depth at the top of the liquefiable interlayer, and H is the thickness of the interlayer. The distribution of the soil layers in each case is shown in Table 1.  In this study, seven different site conditions were constructed, and three sets of numerical models were set up for comparison and analysis according to different variables. The influence of the location of the liquefiable interlayer on the dynamic response of the structure was investigated in cases a, b, and c; the influence of the relative density of the interlayer was investigated in cases b, d, and e; and the influence of the thickness of the interlayer was investigated in cases b, f, and g. The water level is assumed to be at the ground surface in each model, and the liquefiable interlayer is a sandy soil, while the rest of the site soil is non-liquefiable clay. A simplified model of the soil-station structure and the arrangement of the monitoring points are shown in Figure 2. The variable D is the burial depth at the top of the liquefiable interlayer, and H is the thickness of the interlayer. The distribution of the soil layers in each case is shown in Table 1.

Numerical Model Establishment and Parameter Setting
The Finn model was proposed by Martin et al. [19] to solve the problem of volume strain and the pressure change pattern of soil under cyclic loading based on experiments. This model introduces the pore pressure rise model into the Mohr-Coulomb model, which can simulate the process of accumulation and development of super-pore water pressure in soil under dynamic loading until the point of liquefaction. Byrne [20] summarized and analyzed the test data of Martin et al. to simplify the Finn model, and the derived strain increment calculation equation is shown in formula (1):

Numerical Model Establishment and Parameter Setting
The Finn model was proposed by Martin et al. [19] to solve the problem of volume strain and the pressure change pattern of soil under cyclic loading based on experiments. This model introduces the pore pressure rise model into the Mohr-Coulomb model, which can simulate the process of accumulation and development of super-pore water pressure in soil under dynamic loading until the point of liquefaction. Byrne [20] summarized and analyzed the test data of Martin et al. to simplify the Finn model, and the derived strain increment calculation equation is shown in Formula (1): where γ is the shear strain, ε vd is the volume strain, and ∆ε vd is the volume strain increment. In Formula (1), C 1 and C 2 are the constant related to the relative density of the sandy soil and the corrected slamming number, which is calculated empirically using the following formula:

60
(2) Many scholars [21][22][23] have chosen the Finn model to describe the liquefaction behavior of soil and obtained satisfactory data results, which shows that the Finn model is an appropriate tool for investigating the liquefaction principal of soil in this paper.
In order to improve computational accuracy in numerical analysis, interface elements are usually added between the soil and the structure in order to more realistically simulate the nonlinear contact phenomena between the structure and the soil. Empirically, the normal stiffness k n and shear stiffness k s , the key parameters of the contact surface interface, can be taken as 10 times the equivalent stiffness of the surrounding "hardest" adjacent area [24], as shown in Formula (5). The internal friction angle can be approximated as roughly 0.5 times that of the surrounding soil layer: where K is the bulk modulus, G is the shear modulus, and ∆z min is the minimum size of the connection area in the normal direction of the contact surface. Based on the finite difference software FLAC 3D , a numerical model of the interaction between the soil layer and the station structure was established. In order to eliminate the influence of the model boundary effect as much as possible, the calculation range of the foundation is taken as 7 times the length of the station structure, being 150 m in the horizontal direction with a 50 m depth below the ground surface in the vertical direction. The elastic constitutive model is used for concrete, and the Mohr-Coulomb strength criterion is followed for site soils. During the seepage analysis, the soil elements were set up as isotropic models, and the station structure elements were set up as impermeable models. The Byrne-simplified Finn pore pressure model is used to simulate the liquefaction behavior of liquifiable interbedded site soils under seismic action due to the rise in pore water pressure. The physical and mechanical parameters of each stratum soil are given in Table 2, with reference to our engineering experience.
The station is a cast-in-place reinforced concrete structure with a concrete grade of C40, and solid elements are used to model it. To improve the computational efficiency, the numerical model is simplified to a two-dimensional plane strain model considering the interaction between the soil and structure, and the continuity of the central column is considered in light of the equivalent discount of stiffness, which is equated to a continuous longitudinal wall with a thickness of 0.6 m. The interface element is set to simulate the interaction between the soil and the underground station structure. The specific physical and mechanical parameters of the station structure and the interface element parameters are shown in Table 2. The mesh division of the numerical model of the soil-subway station structure for case b is shown in Figure 3

Material Damping Selection
Considering the calculation accuracy and time cost, we used local damping to reflect the nonlinear characteristics of soil materials and the energy loss of the soil-station structure interaction model in the dynamic analysis process. According to our engineering experience and previous research [25,26], the critical damping ratio of rock and soil materials was taken as 10%, and the critical damping ratio of reinforced concrete materials was taken as 5%.
For the sandy soil interlayer, the hysteresis effect due to soil liquefaction needs to be considered further in the dynamic calculation process. Since hysteresis damping can be employed simultaneously with other types of damping and does not affect the calculation time step, hysteresis damping was further applied to the sandy soil interlayer to characterize its hysteresis properties. Among the various models, the Hardin-Drnevich model requires only one parameter for the reference strain, and the fitted equation for the nor- D r : relative density of sandy soil; G: shear modulus; K: bulk modulus; ϕ: friction angle of the soil; ρ d : dry density; n: porosity; k: permeability; c: cohesion; C 1 and C 2 : liquefaction parameters; ρ:density; k s and k n : shear and normal interface stiffness; δ: friction angle of the interface surface.

Material Damping Selection
Considering the calculation accuracy and time cost, we used local damping to reflect the nonlinear characteristics of soil materials and the energy loss of the soil-station structure interaction model in the dynamic analysis process. According to our engineering experience and previous research [25,26], the critical damping ratio of rock and soil materials was taken as 10%, and the critical damping ratio of reinforced concrete materials was taken as 5%.
For the sandy soil interlayer, the hysteresis effect due to soil liquefaction needs to be considered further in the dynamic calculation process. Since hysteresis damping can be employed simultaneously with other types of damping and does not affect the calculation time step, hysteresis damping was further applied to the sandy soil interlayer to characterize its hysteresis properties. Among the various models, the Hardin-Drnevich model requires only one parameter for the reference strain, and the fitted equation for the normalized cut-line modulus (M s ) is shown in Formula (6): For sandy soil materials, the reference strain γ ref is generally taken as 0.06. The modulus decay curves derived from the hysteresis damping model selected for this paper were compared with the experimental data of Seed et al. [27], and the results are shown in Figure 4. It can be seen that the modulus decay curves fitted using the Hardin-Drnevich model are closer to the measured data. were compared with the experimental data of Seed et al. [27], and the results are shown in Figure 4. It can be seen that the modulus decay curves fitted using the Hardin-Drnevich model are closer to the measured data.

Boundary Conditions and Ground Vibration Input
In the process of static analysis, fixed displacement constraint is used at the bottom of the model and in the longitudinal direction, the boundary on both sides is set as a normal displacement constraint, and the ground surface is a free boundary [28]. The initial stress field is obtained after static equilibrium. Then, a natural water level condition is imposed on the ground surface, the bottom and both sides of the model are set as impermeable boundaries, and fluid parameters are established to obtain a constant seepage field. When performing the static analysis, the effects of construction and other factors are ignored, and the ground loads on the underground structure are not considered.
As the bottom element of the model in this numerical analysis model is soft clay, the acceleration time course and velocity time course cannot be applied directly. It is necessary to convert the acceleration and velocity into a stress timescale and then apply them to the bottom of the model. The specific conversion equation [24] is shown in formula (7): where ν s represents the shear components of the velocity at the boundary, ρ is the mass density, and s C represents the s-wave velocities.
Therefore, in the process of dynamic analysis, the existing static constraints at the bottom of the model are first removed, and static boundary conditions are applied at the bottom of the model. Then, the seismic load is applied, and the free field boundary is applied on both sides of the model to reduce the reflection of seismic waves at the model boundary [29,30].
The seismic wave used in this study was the Kobe wave, which is a typical seismic wave measured in structural seismic studies, with an original peak acceleration of 0.85 g in the north-south direction and a duration of the strong seismic part of approximately 7 s. In order to fully reflect the characteristics of site liquefaction, the original peak acceleration of ground shaking is amplitude modulated. The peak acceleration of ground shaking is 0.2 g after amplitude modulation, and the input ground shaking holding time is 20 s. The seismic wave acceleration time curve and Fourier spectrum are shown in Figure 5. To ensure the computational speed and accuracy of the numerical simulation, the amplitude-modulated seismic wave is passed through the fast Fourier transform (FFT) to obtain the input ground shaking power spectrum, the high-frequency components larger than 5

Boundary Conditions and Ground Vibration Input
In the process of static analysis, fixed displacement constraint is used at the bottom of the model and in the longitudinal direction, the boundary on both sides is set as a normal displacement constraint, and the ground surface is a free boundary [28]. The initial stress field is obtained after static equilibrium. Then, a natural water level condition is imposed on the ground surface, the bottom and both sides of the model are set as impermeable boundaries, and fluid parameters are established to obtain a constant seepage field. When performing the static analysis, the effects of construction and other factors are ignored, and the ground loads on the underground structure are not considered.
As the bottom element of the model in this numerical analysis model is soft clay, the acceleration time course and velocity time course cannot be applied directly. It is necessary to convert the acceleration and velocity into a stress timescale and then apply them to the bottom of the model. The specific conversion equation [24] is shown in Formula (7): where ν s represents the shear components of the velocity at the boundary, ρ is the mass density, and C s represents the s-wave velocities. Therefore, in the process of dynamic analysis, the existing static constraints at the bottom of the model are first removed, and static boundary conditions are applied at the bottom of the model. Then, the seismic load is applied, and the free field boundary is applied on both sides of the model to reduce the reflection of seismic waves at the model boundary [29,30].
The seismic wave used in this study was the Kobe wave, which is a typical seismic wave measured in structural seismic studies, with an original peak acceleration of 0.85 g in the north-south direction and a duration of the strong seismic part of approximately 7 s. In order to fully reflect the characteristics of site liquefaction, the original peak acceleration of ground shaking is amplitude modulated. The peak acceleration of ground shaking is 0.2 g after amplitude modulation, and the input ground shaking holding time is 20 s. The seismic wave acceleration time curve and Fourier spectrum are shown in Figure 5. To ensure the computational speed and accuracy of the numerical simulation, the amplitudemodulated seismic wave is passed through the fast Fourier transform (FFT) to obtain the input ground shaking power spectrum, the high-frequency components larger than 5 Hz are filtered out through filtering, and then the filtered seismic wave is baseline-corrected.

Liquefaction Distribution Characteristics
In the numerical calculation process, the excess pore pressure ratio (EPPR) r u is used to represent the liquefaction degree of the soil [24], which is defined as shown in formula (8): where σ ' m is the average effective stress of the element during the dynamic calculation, and σ ' m 0 is the average effective stress of the element before the dynamic calculation. In order to reflect the distribution characteristics of liquefaction zones of liquefiable interlayer sites under different influencing factors, the EPPR clouds of the sandy interlayer after the end of seismic action in different cases are given in Figure 6. When the EPPR is greater than 0.7, the site is considered to be initially liquefied, and when the EPPR is equal to 1, the site is considered to be completely liquefied.
From Figure 6, it can be seen that the liquefaction range of the liquefiable interlayer in each case is basically symmetrically distributed, and the soil near the side wall of the station structure does not liquefy. The reason for this is that the stiffness difference between the station structure and the liquefiable soil layer is too large, and the excess pore water pressure of the soil near the side wall cannot easily accumulate and develop. In case a, when the liquefiable interlayer is located in the upper part of the structure, partial liquefaction occurs in the middle part of the shallow sandy soil on the surface. In case d, when the relative density of the sand interlayer is 30%, almost all the soil in the liquefiable interlayer is liquefied, except for that near the side wall of the structure, and no liquefaction occurs in the liquefiable interlayer in case e, which has the highest relative density. In both case c and case g, liquefaction occurs at the bottom of the structure because the average effective stress at the bottom of the structure is much lower than that at other locations at the same depth due to the presence of the station structure, and liquefaction is more likely to be triggered by the seismic action. In both case b and case f, liquefaction only occurs at the two ends of the negative level of the structure, away from the side walls. In

Liquefaction Distribution Characteristics
In the numerical calculation process, the excess pore pressure ratio (EPPR) r u is used to represent the liquefaction degree of the soil [24], which is defined as shown in Formula (8): where σ m is the average effective stress of the element during the dynamic calculation, and σ m0 is the average effective stress of the element before the dynamic calculation.
In order to reflect the distribution characteristics of liquefaction zones of liquefiable interlayer sites under different influencing factors, the EPPR clouds of the sandy interlayer after the end of seismic action in different cases are given in Figure 6. When the EPPR is greater than 0.7, the site is considered to be initially liquefied, and when the EPPR is equal to 1, the site is considered to be completely liquefied.
From Figure 6, it can be seen that the liquefaction range of the liquefiable interlayer in each case is basically symmetrically distributed, and the soil near the side wall of the station structure does not liquefy. The reason for this is that the stiffness difference between the station structure and the liquefiable soil layer is too large, and the excess pore water pressure of the soil near the side wall cannot easily accumulate and develop. In case a, when the liquefiable interlayer is located in the upper part of the structure, partial liquefaction occurs in the middle part of the shallow sandy soil on the surface. In case d, when the relative density of the sand interlayer is 30%, almost all the soil in the liquefiable interlayer is liquefied, except for that near the side wall of the structure, and no liquefaction occurs in the liquefiable interlayer in case e, which has the highest relative density. In both case c and case g, liquefaction occurs at the bottom of the structure because the average effective stress at the bottom of the structure is much lower than that at other locations at the same depth due to the presence of the station structure, and liquefaction is more likely to be triggered by the seismic action. In both case b and case f, liquefaction only occurs at the two ends of the negative level of the structure, away from the side walls. In summary, the liquefaction of the liquefiable layer is most significant in cases d and g, which will have an important impact on the station structure and should be taken into account in the analysis. summary, the liquefaction of the liquefiable layer is most significant in cases d and g, which will have an important impact on the station structure and should be taken into account in the analysis.

Influence of the Relative Position of the Liquifiable Interlayer
Element A, in the middle of the top plate of the structure, was monitored, and the acceleration time course in each case was obtained as shown in Figure 7. It can be seen Appl. Sci. 2023, 13, 9210 9 of 17 from the figure that the peak acceleration of the top plate of the station is the greatest when the liquifiable interlayer is located on both sides of the station (case b). When the liquifiable interlayer is located at the bottom of the station structure (case c), the peak acceleration of the structure is the lowest, because the soil under the structure is liquefied under the seismic load, and the soil loses its shear characteristics and hinders the transmission of shear waves. Accordingly, the disturbance of the upper metro station is greatly reduced. This indicates that when the liquifiable interlayer is distributed in the lower area of the station, it has a certain seismic isolation effect on the station structure.

Influence of the Relative Position of the Liquifiable Interlayer
Element A, in the middle of the top plate of the structure, was monitored, and the acceleration time course in each case was obtained as shown in Figure 7. It can be seen from the figure that the peak acceleration of the top plate of the station is the greatest when the liquifiable interlayer is located on both sides of the station (case b). When the liquifiable interlayer is located at the bottom of the station structure (case c), the peak acceleration of the structure is the lowest, because the soil under the structure is liquefied under the seismic load, and the soil loses its shear characteristics and hinders the transmission of shear waves. Accordingly, the disturbance of the upper metro station is greatly reduced. This indicates that when the liquifiable interlayer is distributed in the lower area of the station, it has a certain seismic isolation effect on the station structure. The time course curves of the relative horizontal displacement of the top and bottom of the structure when the liquifiable interlayer is located in different spatial positions are given in Figure 8. As can be seen in the figure, the peak relative displacement at the top and bottom of the structure is greatly influenced by the distribution of the liquifiable interlayer location. The peak relative displacement of the structure is the greatest in case b, where the liquifiable interlayer crosses the side wall of the station structure, and the maximum relative displacement is 28.08 mm. Meanwhile, in case c, when the liquifiable interlayer is located at the bottom of the structure, the peak relative displacement is only 12.12 mm, which is reduced by 56.8% compared with case b. From the point of view of structural seismic resistance, when the liquifiable interlayer is located in the lateral foundation of the station structure, the lateral deformation of the structure is the greatest, and it is more likely to be damaged under the earthquake action.
The underground structure is constrained by the surrounding soil layer and mainly undergoes shear deformation under the action of the earthquake. To further investigate the influence of the degree of liquefaction of the interlayer on the structural forces, stress monitoring was conducted, focusing on the key parts of the wall-column structure, and the maximum shear stress at each monitoring point of the structure was obtained as shown in Figure 9. After comparison, it can be seen that the shear stresses at various parts of the structure are higher when the liquifiable interlayer crosses the side wall of the structure (case b), and the shear stresses at the W2 section of the side wall are especially prominent. When the liquifiable interlayer is located at the bottom of the structure (case c), the magnitude of the shear stresses at all parts of the structure is reduced, which also proves that the interlayer has a certain seismic isolation effect when it is located at the bottom of the structure from the perspective of force. The time course curves of the relative horizontal displacement of the top and bottom of the structure when the liquifiable interlayer is located in different spatial positions are given in Figure 8. As can be seen in the figure, the peak relative displacement at the top and bottom of the structure is greatly influenced by the distribution of the liquifiable interlayer location. The peak relative displacement of the structure is the greatest in case b, where the liquifiable interlayer crosses the side wall of the station structure, and the maximum relative displacement is 28.08 mm. Meanwhile, in case c, when the liquifiable interlayer is located at the bottom of the structure, the peak relative displacement is only 12.12 mm, which is reduced by 56.8% compared with case b. From the point of view of structural seismic resistance, when the liquifiable interlayer is located in the lateral foundation of the station structure, the lateral deformation of the structure is the greatest, and it is more likely to be damaged under the earthquake action.
The underground structure is constrained by the surrounding soil layer and mainly undergoes shear deformation under the action of the earthquake. To further investigate the influence of the degree of liquefaction of the interlayer on the structural forces, stress monitoring was conducted, focusing on the key parts of the wall-column structure, and the maximum shear stress at each monitoring point of the structure was obtained as shown in Figure 9. After comparison, it can be seen that the shear stresses at various parts of the structure are higher when the liquifiable interlayer crosses the side wall of the structure (case b), and the shear stresses at the W2 section of the side wall are especially prominent. When the liquifiable interlayer is located at the bottom of the structure (case c), the magnitude of the shear stresses at all parts of the structure is reduced, which also proves that the interlayer has a certain seismic isolation effect when it is located at the bottom of the structure from the perspective of force. Appl. Sci. 2023, 13, x FOR PEER REVIEW 10 of 17

Influence of the Relative Density of the Liquifiable Interlayer
Relative density has an important influence in triggering sand liquefaction. In this section, the most unfavorable case (case b) of the liquifiable interlayer crossing the station structure analyzed in the previous section is used as the base case, and the relative position and thickness of the liquifiable interlayer are kept constant to further investigate the effect of the variation in the relative density of the sand interlayer on the dynamic response of the station structure. Figure 10 shows the acceleration amplification effect of the station structure plate and column. On the whole, the acceleration response of the station structure at different heights differs under the seismic effect, and the peak acceleration of the plate structure is much lower than the peak acceleration of the column structure. The lower the relative density of the sand interlayer is, the stronger the acceleration amplification effect inside the structure will be. The acceleration response at the top and bottom plates of the structure is less affected by the relative density, while the acceleration of the middle plate and the middle column structure is strongly affected by the change in the relative density of the sand interlayer. It is worth noting that the maximum acceleration amplification coefficient of the column structure on the upper level of the station reaches 6.45 in case d. The reason for this is that the great liquefaction of the sandwich soil leads to the weakening of

Influence of the Relative Density of the Liquifiable Interlayer
Relative density has an important influence in triggering sand liquefaction. In this section, the most unfavorable case (case b) of the liquifiable interlayer crossing the station structure analyzed in the previous section is used as the base case, and the relative position and thickness of the liquifiable interlayer are kept constant to further investigate the effect of the variation in the relative density of the sand interlayer on the dynamic response of the station structure. Figure 10 shows the acceleration amplification effect of the station structure plate and column. On the whole, the acceleration response of the station structure at different heights differs under the seismic effect, and the peak acceleration of the plate structure is much lower than the peak acceleration of the column structure. The lower the relative density of the sand interlayer is, the stronger the acceleration amplification effect inside the structure will be. The acceleration response at the top and bottom plates of the structure is less affected by the relative density, while the acceleration of the middle plate and the middle column structure is strongly affected by the change in the relative density of the sand interlayer. It is worth noting that the maximum acceleration amplification coefficient of the column structure on the upper level of the station reaches 6.45 in case d. The reason for this is that the great liquefaction of the sandwich soil leads to the weakening of

Influence of the Relative Density of the Liquifiable Interlayer
Relative density has an important influence in triggering sand liquefaction. In this section, the most unfavorable case (case b) of the liquifiable interlayer crossing the station structure analyzed in the previous section is used as the base case, and the relative position and thickness of the liquifiable interlayer are kept constant to further investigate the effect of the variation in the relative density of the sand interlayer on the dynamic response of the station structure. Figure 10 shows the acceleration amplification effect of the station structure plate and column. On the whole, the acceleration response of the station structure at different heights differs under the seismic effect, and the peak acceleration of the plate structure is much lower than the peak acceleration of the column structure. The lower the relative density of the sand interlayer is, the stronger the acceleration amplification effect inside the structure will be. The acceleration response at the top and bottom plates of the structure is less affected by the relative density, while the acceleration of the middle plate and the middle column structure is strongly affected by the change in the relative density of the sand interlayer. It is worth noting that the maximum acceleration amplification coefficient of the column structure on the upper level of the station reaches 6.45 in case d. The reason for this is that the great liquefaction of the sandwich soil leads to the weakening of the restraint effect of the soil on the structure, which results in greater differential acceleration inside the structure [31,32]. the restraint effect of the soil on the structure, which results in greater differential acceleration inside the structure [31,32]. Figure 10. Acceleration amplification effect of the station plate and column structure. Figure 11 shows the maximum horizontal relative displacement curves of the station side walls under different relative densities of the interlayer. From the figure, it can be seen that the horizontal relative displacement of the structure, with a 70% relative density of the interlayer, is the lowest in case e. As the relative density of the sand interlayer decreases, the degree of liquefaction gradually increases, and the horizontal deformation of the station structure side wall also gradually increases. The maximum relative displacement between the top and bottom plates is 33.4 mm in case d, which increases by 29.4% compared with case e. Figure 11. Influence of the relative density of the liquifiable interlayer on the horizontal relative displacement of the side wall. Figure 12 shows the maximum shear stress at each measured point of the station structure under different relative densities of the liquifiable interlayer. From the figure, it can be seen that the maximum shear stress magnitude at each key section of the wall and column decreases with the increase in the relative density of the sand interlayer, but the degree of reduction varies between different section locations. The shear stress at the  Figure 11 shows the maximum horizontal relative displacement curves of the station side walls under different relative densities of the interlayer. From the figure, it can be seen that the horizontal relative displacement of the structure, with a 70% relative density of the interlayer, is the lowest in case e. As the relative density of the sand interlayer decreases, the degree of liquefaction gradually increases, and the horizontal deformation of the station structure side wall also gradually increases. The maximum relative displacement between the top and bottom plates is 33.4 mm in case d, which increases by 29.4% compared with case e. the restraint effect of the soil on the structure, which results in greater differential acceleration inside the structure [31,32]. Figure 10. Acceleration amplification effect of the station plate and column structure. Figure 11 shows the maximum horizontal relative displacement curves of the station side walls under different relative densities of the interlayer. From the figure, it can be seen that the horizontal relative displacement of the structure, with a 70% relative density of the interlayer, is the lowest in case e. As the relative density of the sand interlayer decreases, the degree of liquefaction gradually increases, and the horizontal deformation of the station structure side wall also gradually increases. The maximum relative displacement between the top and bottom plates is 33.4 mm in case d, which increases by 29.4% compared with case e. Figure 11. Influence of the relative density of the liquifiable interlayer on the horizontal relative displacement of the side wall. Figure 12 shows the maximum shear stress at each measured point of the station structure under different relative densities of the liquifiable interlayer. From the figure, it can be seen that the maximum shear stress magnitude at each key section of the wall and column decreases with the increase in the relative density of the sand interlayer, but the degree of reduction varies between different section locations. The shear stress at the  Figure 11. Influence of the relative density of the liquifiable interlayer on the horizontal relative displacement of the side wall. Figure 12 shows the maximum shear stress at each measured point of the station structure under different relative densities of the liquifiable interlayer. From the figure, it can be seen that the maximum shear stress magnitude at each key section of the wall and column decreases with the increase in the relative density of the sand interlayer, but the degree of reduction varies between different section locations. The shear stress at the bottom of the station structure is less affected by the relative density. The shear stress at the connection between the side wall and the middle plate decreases most significantly with the increase in the relative density and is most affected by the degree of liquefaction of the sand interlayer, which should be given more consideration in seismic design.
Appl. Sci. 2023, 13, x FOR PEER REVIEW 12 of 17 bottom of the station structure is less affected by the relative density. The shear stress at the connection between the side wall and the middle plate decreases most significantly with the increase in the relative density and is most affected by the degree of liquefaction of the sand interlayer, which should be given more consideration in seismic design. Figure 12. Influence of the relative density of the liquifiable interlayer on the shear stress at each section. Figure 13 shows the time course curve of the soil and water compressive stress in the middle of the side-wall-monitoring element S. It can be seen that the development trend of the soil and water compressive stress is similar in each case, always progressing through the three stages of "smooth development, fluctuating rise and stabilization". The lateral wall soil and water compressive stresses show a decreasing trend with the increase in the relative density of the sand and soil interlayer, and the lower the relative density is, the greater the peak soil and water compressive stresses in the lateral wall are. Because the relative density of the sandwich soil decreases, the degree of liquefaction increases, the development of super-pore water pressure is more significant, and the soil and water compressive stress on the side wall increases, which further increases the influence of horizontal shear on the side wall.   Figure 13 shows the time course curve of the soil and water compressive stress in the middle of the side-wall-monitoring element S. It can be seen that the development trend of the soil and water compressive stress is similar in each case, always progressing through the three stages of "smooth development, fluctuating rise and stabilization". The lateral wall soil and water compressive stresses show a decreasing trend with the increase in the relative density of the sand and soil interlayer, and the lower the relative density is, the greater the peak soil and water compressive stresses in the lateral wall are. Because the relative density of the sandwich soil decreases, the degree of liquefaction increases, the development of super-pore water pressure is more significant, and the soil and water compressive stress on the side wall increases, which further increases the influence of horizontal shear on the side wall. bottom of the station structure is less affected by the relative density. The shear stress at the connection between the side wall and the middle plate decreases most significantly with the increase in the relative density and is most affected by the degree of liquefaction of the sand interlayer, which should be given more consideration in seismic design. Figure 12. Influence of the relative density of the liquifiable interlayer on the shear stress at each section. Figure 13 shows the time course curve of the soil and water compressive stress in the middle of the side-wall-monitoring element S. It can be seen that the development trend of the soil and water compressive stress is similar in each case, always progressing through the three stages of "smooth development, fluctuating rise and stabilization". The lateral wall soil and water compressive stresses show a decreasing trend with the increase in the relative density of the sand and soil interlayer, and the lower the relative density is, the greater the peak soil and water compressive stresses in the lateral wall are. Because the relative density of the sandwich soil decreases, the degree of liquefaction increases, the development of super-pore water pressure is more significant, and the soil and water compressive stress on the side wall increases, which further increases the influence of horizontal shear on the side wall.

Influence of the Layer Thickness of the Liquifiable Interlayer
The analysis in this section, again, takes case b as the benchmark case, keeps the relative position and relative density of the liquifiable interlayer unchanged, gradually increases the thickness of the sand layer, and further explores the influence of the thickness change in the liquifiable interlayer on the dynamic response of the station structure. Figure 14 shows the maximum horizontal displacement curve of the side wall of the station under different layer thicknesses of the liquifiable interlayer. It can be seen from the figure that when the whole liquifiable interlayer is only distributed in the lateral foundation of the structure (case b and case f), as the thickness of the interlayer increases, the horizontal relative displacement amplitude of the side walls gradually increases, and the lateral deformation of the structure also increases. When the thickness of the liquifiable layer is 12 m (case f), the relative displacement between the top and bottom of the structure is the greatest, reaching 46.6 mm, which is 1.6 times that of the site where the interlayer thickness is 5 m. In case g, when the thickness of the liquifiable layer is 20 m, the station structure is completely placed in the liquefiable interlayer, and the relative displacement of the side wall of the structure is the lowest, only 26.5 mm. Under earthquake action, the water and soil pressure produced by the liquefiable sand interlayer on the structure is greater than that of the non-liquefiable clay layer. As the thickness of the liquefiable layer increases, the contact area between the side wall of the station and the liquefiable interlayer increases, further increasing the lateral deformation of the structure; When the structure is completely placed in the liquefiable layer, the water and soil pressure difference on the structure is greatly reduced, so the deformation of the structure is minimal.

Influence of the Layer Thickness of the Liquifiable Interlayer
The analysis in this section, again, takes case b as the benchmark case, keeps the relative position and relative density of the liquifiable interlayer unchanged, gradually increases the thickness of the sand layer, and further explores the influence of the thickness change in the liquifiable interlayer on the dynamic response of the station structure. Figure 14 shows the maximum horizontal displacement curve of the side wall of the station under different layer thicknesses of the liquifiable interlayer. It can be seen from the figure that when the whole liquifiable interlayer is only distributed in the lateral foundation of the structure (case b and case f), as the thickness of the interlayer increases, the horizontal relative displacement amplitude of the side walls gradually increases, and the lateral deformation of the structure also increases. When the thickness of the liquifiable layer is 12 m (case f), the relative displacement between the top and bottom of the structure is the greatest, reaching 46.6 mm, which is 1.6 times that of the site where the interlayer thickness is 5 m. In case g, when the thickness of the liquifiable layer is 20 m, the station structure is completely placed in the liquefiable interlayer, and the relative displacement of the side wall of the structure is the lowest, only 26.5 mm. Under earthquake action, the water and soil pressure produced by the liquefiable sand interlayer on the structure is greater than that of the non-liquefiable clay layer. As the thickness of the liquefiable layer increases, the contact area between the side wall of the station and the liquefiable interlayer increases, further increasing the lateral deformation of the structure; When the structure is completely placed in the liquefiable layer, the water and soil pressure difference on the structure is greatly reduced, so the deformation of the structure is minimal. The shear stresses in each key part of the structure under different thicknesses of liquefaction interlayer are shown in Figure 15. From the figure, it can be seen that when all the liquifiable interlayer is only located in the lateral foundation of the structure (case b and case f), the influence law of the liquifiable layer thickness on each part of the station structure is not consistent. The variation of liquifiable interlayer thickness does not have much effect on the magnitude of shear stress in each section of the central column, the shear stress at the intersection of the structural side wall and the top and bottom plate gradually increases with the increase in the interlayer thickness, but the magnitude of the shear stress at the connection between the side wall and the central plate gradually decreases. When the thickness of the liquifiable interlayer is as great as 20 m (case g), the amplitude of the shear stress in all parts of the station wall and column structure is the The shear stresses in each key part of the structure under different thicknesses of liquefaction interlayer are shown in Figure 15. From the figure, it can be seen that when all the liquifiable interlayer is only located in the lateral foundation of the structure (case b and case f), the influence law of the liquifiable layer thickness on each part of the station structure is not consistent. The variation of liquifiable interlayer thickness does not have much effect on the magnitude of shear stress in each section of the central column, the shear stress at the intersection of the structural side wall and the top and bottom plate gradually increases with the increase in the interlayer thickness, but the magnitude of the shear stress at the connection between the side wall and the central plate gradually decreases. When the thickness of the liquifiable interlayer is as great as 20 m (case g), the amplitude of the shear stress in all parts of the station wall and column structure is the lowest, being much lower than in other cases. The reason for this phenomenon is that the soil in the lower part of the structure liquefies first, and the soil near the sides of the side walls flows downward, which slows down the effect on the horizontal shear of the structure. lowest, being much lower than in other cases. The reason for this phenomenon is that the soil in the lower part of the structure liquefies first, and the soil near the sides of the side walls flows downward, which slows down the effect on the horizontal shear of the structure. The research of Xia et al. [33] showed that the thickness of the liquifiable soil layer is one of the key factors affecting the uplift of underground structures. Firstly, the developmental trend of the vertical displacement of the structure is analyzed. Taking case b as an example, it can be seen from Figure 16 that the structure gradually starts to float within 0~4 s, and the floating amplitude is almost negligible. With the increase in seismic intensity, the vertical displacement of each monitoring point climbs rapidly for around 4 s and then reaches the peak and fluctuates; at the end of the seismic action, the vertical displacement tends to become stable and remains basically unchanged. To further investigate the uplift response of the station structure under different liquifiable interlayer thicknesses, the final uplift displacement values of the structure under different layer thicknesses of the liquifiable interlayer are given in Figure 17. It can be seen that the floating displacement of the station structure increases with the increase in the The research of Xia et al. [33] showed that the thickness of the liquifiable soil layer is one of the key factors affecting the uplift of underground structures. Firstly, the developmental trend of the vertical displacement of the structure is analyzed. Taking case b as an example, it can be seen from Figure 16 that the structure gradually starts to float within 0~4 s, and the floating amplitude is almost negligible. With the increase in seismic intensity, the vertical displacement of each monitoring point climbs rapidly for around 4 s and then reaches the peak and fluctuates; at the end of the seismic action, the vertical displacement tends to become stable and remains basically unchanged. lowest, being much lower than in other cases. The reason for this phenomenon is that the soil in the lower part of the structure liquefies first, and the soil near the sides of the side walls flows downward, which slows down the effect on the horizontal shear of the structure. The research of Xia et al. [33] showed that the thickness of the liquifiable soil layer is one of the key factors affecting the uplift of underground structures. Firstly, the developmental trend of the vertical displacement of the structure is analyzed. Taking case b as an example, it can be seen from Figure 16 that the structure gradually starts to float within 0~4 s, and the floating amplitude is almost negligible. With the increase in seismic intensity, the vertical displacement of each monitoring point climbs rapidly for around 4 s and then reaches the peak and fluctuates; at the end of the seismic action, the vertical displacement tends to become stable and remains basically unchanged. To further investigate the uplift response of the station structure under different liquifiable interlayer thicknesses, the final uplift displacement values of the structure under different layer thicknesses of the liquifiable interlayer are given in Figure 17. It can be seen that the floating displacement of the station structure increases with the increase in the To further investigate the uplift response of the station structure under different liquifiable interlayer thicknesses, the final uplift displacement values of the structure under different layer thicknesses of the liquifiable interlayer are given in Figure 17. It can be seen that the floating displacement of the station structure increases with the increase in the thickness of the liquifiable interlayer, and there is a clear difference in the floating phenomenon. The specific performance is such that the floating displacement on the left side is lower than that on the right. When the station structure is wholly placed in the liquifiable layer (case g), the degree of floating is the greatest. Although the shearing effect on the structure in the horizontal direction is weakened at this time, the adverse effect on the structure due to the floating effect still needs to be considered in the seismic design.
Appl. Sci. 2023, 13, x FOR PEER REVIEW 15 of 17 thickness of the liquifiable interlayer, and there is a clear difference in the floating phenomenon. The specific performance is such that the floating displacement on the left side is lower than that on the right. When the station structure is wholly placed in the liquifiable layer (case g), the degree of floating is the greatest. Although the shearing effect on the structure in the horizontal direction is weakened at this time, the adverse effect on the structure due to the floating effect still needs to be considered in the seismic design.

Conclusions
Based on the finite difference numerical calculation method, the influence law of the liquifiable interlayer on the dynamic response of underground station structures was analyzed, and the influences of interlayer location distribution, relative density, and layer thickness on structural force deformation were explored. The main conclusions are as follows: (1) The liquefaction range of the liquefiable interlayer is basically symmetrically distributed, and the liquefied area is mainly distributed at the bottom of the station structure and on both sides, away from the structure area. When the wall stiffness is great compared with the surrounding soil, the site soil near the structural side wall cannot easily liquefy. (2) The relative location distribution of the liquifiable interlayer has a strong influence on the structural seismic response. When the liquifiable interlayer crosses the middle of the station, the shear stress and lateral deformation of the structure are the greatest, while when the liquifiable interlayer is located in the bottom area of the station structure, the dynamic response of the station structure is greatly reduced, and at this time, the liquifiable interlayer has certain effects of seismic damping and seismic isolation, which are beneficial in supporting the seismic resistance of the station structure to a certain extent. (3) When the liquifiable interlayer is distributed in the middle of the lateral foundation of the station structure, the shear stress of the negative one-story wall and column decreases with increasing relative density, and the structural shear stress of the wall column near the bottom plate is not greatly affected by the relative density. When the relative density of the sand interlayer is low, the internal acceleration response of the structure is strong, and the acceleration amplification effect of the middle column structure is especially prominent. (4) When the whole liquifiable interlayer is only located in the lateral foundation of the structure, with an increase in the thickness of the interlayer, the shearing effect on the

Conclusions
Based on the finite difference numerical calculation method, the influence law of the liquifiable interlayer on the dynamic response of underground station structures was analyzed, and the influences of interlayer location distribution, relative density, and layer thickness on structural force deformation were explored. The main conclusions are as follows: increases to a certain extent, the station structure is wholly placed in the liquifiable interlayer, and although the horizontal disturbance of the structure decreases at this time, the uplift displacement of the structure reaches its maximum.
Author Contributions: Writing-original draft preparation, J.Y.; writing-review and editing, Y.L. and J.Y. All authors have read and agreed to the published version of the manuscript.