As one of the important types of forest resources, mapping forest stock volume (FSV) in larch (Larix decidua) forests holds significant importance for forest resource management, carbon cycle research, and climate change monitoring. However, the accuracy of FSV mapping using common spectral and texture features is often limited due to their failure in fully capturing seasonal changes and growth cycle characteristics of vegetation. Phenological features can effectively provide essential information regarding the growth status of forests. In this study, multi-temporal Sentinel-2 satellite imagery were initially acquired in the Wangyedian Forest Farm in Chifeng City, Inner Mongolia. Subsequently, various phenological features were extracted from time series variables constructed by Gaussian Process Regression (GPR) using Savitzky–Golay filters, stepwise differentiation, and Fourier transform techniques. The alternative features were further refined through Pearson’s correlation coefficient analysis and the forward selection algorithm, resulting in six groups of optimal subsets. Finally, four models including the Random Forest (RF), K-Nearest Neighbors (KNN), Support Vector Machine (SVM), and Multiple Linear Regression (MLR) algorithms were developed to estimate FSV. The results demonstrated that incorporating phenological features significantly enhanced model performance, with the SVM model exhibiting the best performance—achieving an R² value of 0.77 along with an RMSE value of 46.36 m³/hm² and rRMSE value of 22.78%. Compared to models without phenological features, inclusion of these features led to a 0.25 increase in R² value while reducing RMSE by 10.40 m³/hm² and rRMSE by 5%. Overall, integration of phenological feature variables not only improves the accuracy of larch forest FSV mapping but also has potential implications for delaying saturation phenomena.