【Objective】The argillic horizon is a subsurface secondary layer formed by the accumulation of soil clay particles, and its thickness exerts a crucial regulatory effect on soil processes and vegetation growth in Alfisols. Understanding its spatial distribution is critical for effective land management, particularly in agriculturally important regions such as Northeast China. However, there is still limited knowledge of the spatial variability in argillic horizon thickness, and predictive studies on this topic are scarce. Traditional understanding has largely relied on extensive field surveys combined with geostatistical methods, which are often resource-intensive and may not be efficient over large regions. This study aims to develop a robust predictive model to map the spatial distribution of argillic horizon thickness across the three northeastern provinces of China by integrating limited soil profile observations with a rich set of environmental covariates. 【Method】A total of 311 soil profile samples with argillic horizons were collected in Northeast China. These samples incorporated data from recent field surveys and historical soil records. Consistent with the SCORPAN framework, 71 environmental covariates were selected to correspond to relief, climate, organism, and soil factors. Dual feature selection was conducted via Pearson correlation analysis and the Boruta algorithm. The quantile regression forest (QRF) model was then adopted for spatial modeling, cross-validation, and uncertainty estimation. Rigorous evaluation of model performance and uncertainty estimation was conducted through 50 repetitions of 10-fold cross-validation, and accumulated local effects (ALE) plots were generated to interpret the relationship between key predictors and the target variable. 【Result】The average results from 50 iterations showed that the model achieved a coefficient of determination (R2) of 0.32, a root mean square error (RMSE) of 24.34 cm, and a mean absolute error (MAE) of 19.47 cm. This performance is significantly superior to that of most regional and national scale soil thickness prediction studies (R2 = 0.11–0.41). The prediction interval coverage percentage (PICP) was 86.2%, which is close to the predefined 90% prediction interval (PI), indicating high reliability of the uncertainty estimation. Soil and climate factors were generally more influential than organism and relief factors, with soil thickness (ST) identified as the most critical driving factor. The spatial prediction results indicated a distinct decreasing trend in argillic horizon thickness from the southwest to the northeast. The western and southwestern regions of the study area exhibited the thickest argillic horizons (mostly over 80 cm, with some regions ranging from 100 to 125 cm), while the northern, eastern, and southeastern regions had thinner ones (mostly 20–35 cm, with some regions below 20 cm). High prediction uncertainty was concentrated in mountainous and hilly regions with sparse soil survey points. 【Conclusion】This study confirms the feasibility of mapping argillic horizon thickness using a machine learning approach combined with environmental covariates, even in large, complex landscapes with limited soil observations. Future research could focus on integrating proxies for parent material and pedogenic age to enhance model accuracy, as well as exploring the spatial prediction of other argillic horizon properties (e.g., upper boundary and compactness). This study not only addresses the gap in argillic horizon thickness prediction in Northeast China, but also offers valuable insights for optimizing regional land management strategies.