Snow density is one of the essential properties to describe snowpack characteristics. To obtain the spatial variability of snow density and estimate it accurately in different periods of snow season still remain as challenges, particularly in the mountains. This study analyzed the spatial variability of snow density with in-situ measurements in three different periods (i.e. accumulation, stable, melt period) of snow seasons 2017/2018 and 2018/2019 in the middle Tianshan Mountains, China. The performance of multiple linear regression model (MLR) and three machine learning models (i.e. Random Forest (RF), Extreme Gradient Boosting (XGB) and Light Gradient Boosting Machine (LGBM)) to simulate snow density has been evaluated. It was found that the snow density in melt period (0.27 g cm-3) was generally greater than that in stable (0.20 g cm-3) and accumulation period (0.18 g cm-3), and the spatial variability of snow density in melt period was slightly smaller than that in other two periods. The snow density in mountainous areas was generally higher than that in plain or valley areas, and snow density increased significantly (p < 0.05) with elevation in the accumulation and stable periods. Besides elevation, latitude and ground surface temperature also had critical impacts on the spatial variability of snow density in the middle Tianshan Mountains, China. In this work, the machine learning model, especially RF model, performed better than MLR on snow density simulation in three periods. Compared with MLR, the determination coefficients of RF promoted to 0.61, 0.51 and 0.58 from 0.50, 0.1 and 0.52 in accumulation period, stable period and melt period respectively. This study provide a more accurate snow density simulation method for estimating regional snow mass and snow water equivalent, which allows us to achieve a better understanding of regional snow resources.