Information about forest stand species distribution is essential for biodiversity modelling, forest disturbances, fire hazard and drought monitoring, biomass and carbon estimation, detection of non-native and invasive species, as well as for planning forest management strategies. High temporal and spectral resolution remote sensing data from the Sentinel-2 mission enables the derivation of accurate and timely maps of tree species in forests in a cost-efficient way. However, there is still a lack of studies regarding forest stand species mapping for large areas like the Polish Carpathian Mountains (approx. 20,000 km2). In this study, we aimed to develop a workflow to obtain forest stand species maps with machine learning algorithms applied to multi-temporal Sentinel-2 products and environmental data at regional scale. Using variable importance techniques - Variable Importance Using Random Forests (VSURF) and Recursive Feature Elimination (RFE) - we assessed three Sentinel-2 Best Available Pixel composites (April, July and October), eight annual spectral-temporal metrics (STM; mean, minimum, maximum, standard deviation, range, first quartile, third quartile and interquartile range), and four environmental topographic variables (elevation, slope, aspect, distance to water bodies), i.e. 114 variables in total. Following a variable importance assessment, we produced maps of eleven tree species with the use of three Machine Learning algorithms: Random Forest (RF), Support Vector Machines (SVM) and Extreme Gradient Boosting (XGB) on nine different variable subsets, i.e. in total 27 classifications. The results showed that SVM outperformed the other two classifiers - the highest overall accuracy exceeded 85% for SVM classification of all variables (86.9%), and 64 variables (85.6%). Including elevation information improved the accuracies. From the best five classifications we created a final ensemble map (overall accuracy of 86.6%) and a precision map based on the Simpson Index, which indicates where the five models agree. This ensemble approach allowed us to determine that the lowest precision occurred in foothills and basins with lower forest cover, in the areas with lack of good quality imagery, and at the borders of stands with homogenous species composition. On the other hand, the highest precision occurred in regions with homogenous stands with high forest and canopy cover. The study demonstrates the potential of Sentinel-2 imagery and topographic data in mapping forest stand species in large mountainous areas with high accuracy. Furthermore, it demonstrates the usefulness of the ensemble approach, which enables to assess the classification precision.