In the construction process of the artificial ground freezing (AGF), the utilization of the temperature field to determine the freezing time is crucial for the safe construction. While the thermal parameter is the core parameter of the temperature field calculation. How to obtain the joint probability distribution of thermal parameters of frozen soil under limited test data is essential to improve the accuracy of the stochastic temperature field and guide safe construction. In this study, based on the sample of frozen soil in Mengtie railway station of Anhui metro line, the statistical characteristics of the thermal conductivity, volumetric heat capacity, and thermal conductivity at various temperatures were obtained. Then multidimensional Gaussian copula models at different temperatures were constructed using three construction methods to characterize the correlations among the thermal parameters. Additionally, the correlation coefficients under three methods were simulated using Monte Carlo simulation (MCS) and the fitting errors were calculated. Finally, the Sobol indices of thermal parameters were calculated by a simple one-dimensional heat conduction model. The results show that the frozen soil thermal parameters have obvious correlation variability at various temperatures. The Pearson method presents the most favorable fitting capability in the construction of the joint distribution model of frozen soil thermal variables. The error of the simulation results is the smallest when the construction method and correlation coefficient are identical. The Sobol indices calculated by different methods have significant differences, with the Sobol indices using the Kendall method exhibiting a higher sensitivity to the nonlinearity of the parameters.