Multi-model ensemble mean shows accelerating global below-ground warming

The below-ground component of the Earth’s critical zone is crucial to human activities and underpins numerous chemical, physical and biological processes. However, previous studies primarily concentrated on below-ground temperatures (BGT) until depths up to 3 m and periods shorter than 70 years; few studies have globally analyzed the historical spatiotemporal variability of BGT beyond those ranges. The objective of this study was to investigate BGT anomalies (ΔBGT) between depths of 0–42 m during 1850–2100 using model outputs from CMIP6. The results show a three-stage accelerating warming pattern (1850–2014): weak pre-1943 warming (0.02 °C decade⁻¹, depth-average), mid-century stagnation, and post-1984 acceleration (0.33 °C decade⁻¹, depth-average) for depth mean of 0.05∼1.75 m. Future mean warming rises ∼1.7 times from SSP1‑2.6 (2.08 °C) to SSP5‑8.5 (3.45 °C), with maximum of warming mean expanding 2.6 times. Asymmetric BGT extremes drive elevated subsurface heat risk under high emissions. A robust seasonal hierarchy reversal occurs (DJF‑ to JJA‑dominated), with winter BGT most sensitive to radiative forcing. ΔBGT amplifies strongly from 60°N, and enhances in high‑altitude/coastal regions under high emissions. Heterogeneous bottom boundary condition placement (BBCP) is an important structured uncertainty source in multi-model BGT analysis, introducing non-physical sampling artifacts in ensemble-mean vertical profiles. Despite inter‑model heterogeneity, the multi‑model ensemble yields physically consistent depth‑attenuated warming, providing an ensemble-constrained reference for subsurface thermal change investigation. By 2100, low-moderate emission scenarios (e.g., SSP1‑2.6, SSP2‑4.5) will slow BGT warming. This study can provide insightful understanding of the overlooked BGT and inform future model intercomparison projects and ensemble mean analysis.