## 6 Analysis of Model Parameters  
In this section, the control variate method [32; 33] is employed to explore the influences of the main model parameters, such as the external pressure \(F\_0\), influence coefficient \(b_{\mathrm{r}}\), the radial gradient influence coefficient \(\alpha_{\mathrm{f}}\), and the anisotropic frost heave coefficient \(k_{\mathrm{f}}\) on the stress, displacement, and plastic zone radius, and frost-heaving force. The values of the model parameters used for the control variate analysis are listed in Table 1.  
The stress distribution curves of the rock mass and the support structure under different model parameters are shown in Figure 10. Regarding the radial stress \(\sigma_{r}\), except for the slight effect of \(b_{\mathrm{r}}\) on \(\sigma_{r}\) in the plastic zone of the frozen rock mass II-1, the influences of parameters \(F\_0\), \(\alpha_{\mathrm{f}}\), and \(k_{\mathrm{f}}\) on \(\sigma_{r}\) were negligible, and the maximum deviation between different parameters was less than 2.0%. Regarding the circumferential stress \(\sigma_{\theta}\), the impacts of parameters \(F\_0\), \(b_{\mathrm{r}}\), \(\alpha_{\mathrm{f}}\), and \(k_{\mathrm{f}}\) on \(\sigma_{\theta}\) were mainly reflected in support zone I and frozen rock mass zone II, and the effects in unfrozen rock mass III were likewise negligible. As \(F\_0\) increased from 0.0 MPa to 1.0 MPa, \(\sigma_{\theta}\) increased in zones I and II-1 and decreased in zone II-2. The influence of the intensity of \(F\_0\) on \(\sigma_{\theta}\) in zone I was greater than that in zones II-1 and II-2. As \(b_{\mathrm{r}}\) increased from 0.0 to 1.0, \(\sigma_{\theta}\) decreased in zone I and increased in zones II-1 and II-2. The influence of the intensity of \(b_{\mathrm{r}}\) on \(\sigma_{\theta}\) in zones I and II-1 was greater than that in zone II-2, especially in zone I and the elastic-plastic contact surface in zone II. As parameter \(\alpha_{\mathrm{f}}\) increased from 0.0 to 0.4, \(\sigma_{\theta}\) increased in zones I and II-1 and decreased in zone II-2. The influence of the intensity of \(\alpha_{\mathrm{f}}\) on \(\sigma_{\theta}\) in zone I was greater than that in zones II-1 and II-2. As \(k_{\mathrm{f}}\) increased from 1.0 to 3.0, \(\sigma_{\theta}\) increased in zones I and II-1 and decreased in zone II-2. The influence of the intensity of parameters \(k_{\mathrm{f}}\) on \(\sigma_{\theta}\) in zone II-2 was greater than that in zones I and II-1. Figure 10: Stress distribution under different model parameters: (**a**) pressure acting on inner surface of support \(F\_0\); (**b**) influence coefficient of intermediate principal stress \(b_{\mathrm{r}}\); (**c**) radial gradient influence coefficient \(\alpha_{\mathrm{f}}\); (**d**) anisotropic frost heave coefficient \(k_{\mathrm{f}}\). The effects of the model parameters on the plastic zone radius \(R_{\rm P}\) and the frost-heaving force \(F\_1\) are shown in Figure 11. \(R_{\rm P}\) decreased with the increase in \(F\_0\), \(b_{\rm r}\), \(\alpha_{\rm f}\), and \(k_{\rm f}\). A significant negative correlation exists between \(R_{\rm P}\) and \(F\_0\), \(b_{\rm r}\), \(\alpha_{\rm f}\), and \(k_{\rm f}\). For example, when \(b_{\rm r}\) increased from 0.25 to 0.75, \(R_{\rm P}\) decreased from 5.86 m to 5.76 m, and the decreasing amplitude was approximately 1.7%. When \(\alpha_{\rm fr}\) increased from 0.0 to 0.4, \(R_{\rm P}\) decreased from 5.94 m to 5.87 m, and the decrease in amplitude was approximately 1.2%. Moreover, \(F\_1\) increased with the increase in \(F\_0\), \(\alpha_{\rm f}\), and \(k_{\rm f}\), and a positive correlation exists between \(F\_1\) and \(F\_0\), \(\alpha_{\rm f}\), and \(k_{\rm f}\). \(F\_1\) decreased with the increase in \(b_{\rm r}\), and \(F\_1\) is negatively correlated with \(b_{\rm r}\). For example, when \(b_{\rm r}\) increased from 0.25 to 0.75, \(F\_0\) decreased from 1.96 MPa to 1.54 MPa, and the decreasing amplitude was approximately 27.3%.