Ice-rich permafrost thaws in response to rapid Arctic warming, and ground subsidence facilitates the formation of thermokarst lakes. Thermokarst lakes transform the surface energy balance of permafrost, affecting geomorphology, hydrology, ecology, and infrastructure stability, which can further contribute to greenhouse gas emissions. Currently, the spatial distribution of thermokarst lakes at large scales remains a challenging task. Based on multiple high-resolution environmental factors and thermokarst lake inventories, we used machine learning methods to estimate the spatial distributions of present and future thermokarst lake susceptibility (TLS) maps. We also identified key environmental factors of the TLS map. At 1.8 × 10⁶ km², high and very high susceptible regions were estimated to cover about 10.4 % of the region poleward of 60°N, which were mainly distributed in permafrost-dominated lowland regions. At least 23.9 % of the area of TLS maps was projected to disappear under representative concentration pathway scenarios (RCPs), with increased susceptibility levels in northern Canada. The slope was the key conditioning factor for the occurrence of thermokarst lakes in Arctic permafrost regions. Compared with similar studies, the reliability of the TLS map was further evaluated using probability calibration curve and coefficient of variation (CV). Our results provide a means for assessing the spatial distribution of thermokarst lakes at the circum-Arctic scale but also improve the understanding of their dynamics in response to the climate system.