In the context of global warming, understanding the evolving patterns of extreme snowfall is a critical scientific challenge. This study investigates the spatiotemporal trends, atmospheric drivers, and thermodynamic mechanisms of snowfall in China from 1961 to 2020. Precipitation phases were classified using the highly accurate Ding's wet-bulb temperature method. A T-Learner meta-learning approach was used to determine the causal impact of atmospheric circulation, and the Clausius-Clapeyron (CC) relationship was examined to understand temperature dependency. The results indicated that: 1) The Ding's wet-bulb temperature method proves to be a robust tool for classifying historical precipitation phases in China. Its accuracy in snowfall detection (POD) and amount estimation (R) is high nationwide but shows a clear geographical gradient, with performance systematically improving from south to north. 2) China's snowfall exhibited a paradoxical trend, characterized by "increased total snowfall, reduced snowfall days, intensified extreme events, shorter snowfall duration, and longer snow-free periods." The snowfall amount (ST, +0.31 mm/decade) significant increases and snowfall days (SD, -0.29 days/decade) significant decreases. Maximum consecutive snowfall days (MCSD, -0.050 days/ decade) slight decrease and maximum no-snowfall days (MNSD, +1.50 days/decade) slight increase. Some extreme snowfall metrics Maximum daily (Sx1day, +0.16 mm/decade) slight increase multi-day snowfall (Sx3day, +0.35 mm/decade) significant increases. Extreme event intensity 90th percentile of snowfall amount (S90p, +0.45 mm/decade), 95th percentile of snowfall amount (S95p, +0.37 mm/decade) and days exceeding these percentiles (S90pD, +0.057 days/decade), (S95pD, +0.036 days/decade) all showed significant increases. 3) Using a T-Learner causal inference framework with XGBoost, LightGBM, and Random Forest as base models, we investigated the causal impacts on extreme snowfall (S95p). We found that a single circulation index can have opposing effects (promoting vs. suppressing) across different regions. The Pacific-North American (PNA) pattern was identified as the dominant positive driver in most regions, whereas decadal modes like the Pacific Decadal Oscillation (PDO) and Atlantic Multidecadal Oscillation (AMO) typically exert a suppressive influence. This finding was consistently confirmed across all three base models, underscoring the robustness of our results. 4) Across -20 degrees C to 0 degrees C, all subregions exhibit a consistent sub-Clausius-Clapeyron scaling, with the semi-humid zone (SHZ) showing the strongest warming sensitivity (+5.4 %/degrees C) and the semi-arid zone (SAZ) the weakest (+2.7 %/degrees C). Above 0 degrees C, extreme snowfall intensity declines sharply with temperature, exhibit a negative Clausius-Clapeyron scaling relationship. The humid zone (HZ) shows the largest decrease (-20.5 %/degrees C), whereas the SAZ response is minimal (-0.2 %/degrees C).