1. Introduction

Extreme low-temperature events pose significant climate hazards by causing freezing rain, snowstorms, power disruptions [1], crop damage that threatens food security [2], and increased cardiovascular disease risks [3]. Understanding their spatiotemporal patterns and driving mechanisms remains a key focus in atmospheric science. Located in the East Asian monsoon zone, China's climate exhibits exceptional complexity due to interactions among Eurasia-Pacific-Indian Ocean thermal contrasts and the Tibetan Plateau's orographic forcing [4], resulting in frequent extreme weather. Despite rising winter temperatures (0.8–1.2°C) in Northeast and East China under East Asian winter monsoon modulation and global warming [5,6], extreme cold events paradoxically show intensified frequency fluctuations [7]. Notably, their climate sensitivity index (CSI) exceeds 1.7 – higher than for extreme high-temperature events [8] – suggesting amplified impacts of global warming.

Identifying changing patterns is crucial for revealing driving factors. Arctic amplification and atmospheric circulation anomalies synergistically intensify these events, though with significant regional heterogeneity [9]. Northwest China shows marked reductions post-1985 [10], while North China exhibits decadal shifts: increased extreme cold days after 2003 followed by decreases after 2013 [11]. Northeast China displays a triphasic trend (high frequency pre-1980s, decline until early-2000s, resurgence post-2000) [9,12]. This spatial variability aligns with amplified warming gradients across China [13]. Paradoxically, while frequency decreases under global warming, intensity has increased [8]. Both intrinsic atmospheric variability and external diabatic forcing (e.g., sea ice anomalies, SST variations) contribute [14].

Existing studies on China's winter extreme cold events primarily focus on local analyses, lacking comprehensive examination of nationwide spatial modes, interannual-decadal evolution patterns, and circulation-driven mechanisms. Here, we combine high-resolution reanalysis data with statistical-physical hybrid methods to clarify multi-scale spatiotemporal characteristics of extreme low-temperature days across China (1961–2020) and reveal their physical links to atmospheric circulation and SSTA anomalies.

2. Data and methods

2.1. Long time series reanalysis data

This study used observation and reanalysis data, including: CN05.1 daily minimum temperature provided by the National Meteorological Information Center , which is based on the observation data from more than 2,400 ground meteorological stations in China. Using the anomaly approximation method, it is obtained by interpolating and superimposing the climate field and anomaly field separately [15], with a resolution of 0.25°×0.25°; the reanalysis dataset released by the National Center for Environmental Prediction/National Center for Atmospheric Research (NCEP/NCAR ReanalysisⅠ) [16], which includes the global monthly sea level pressure field (SLP), 850 hPa, 500 hPa, and 200 hPa geopotential height field and wind field data with a spatial resolution of 2.5°×2.5°; the global sea surface temperature dataset developed by the Met Office Hadley Centre (HadISST) with a spatial resolution of 1°×1°. The winter referred to in the article is from December to February of the following year, which is DJF. After screening the data, a total of 60 winters were identified from 1961 to 2020.

2.2. Construction of extreme low temperature days dataset

Referring to the definition of extreme low temperature thresholds [17-19], this paper used the CN05.1 daily minimum temperature data, for a given grid point, all the daily minimum temperatures for the 60 years from 1961 to 2020 (21,915 days in total) are sorted in ascending order, and the temperature value at the fifth percentile is taken as the extreme low temperature threshold for that grid point. If the daily minimum temperature on a certain day in that grid point falls below the threshold, it is assumed that an extreme low temperature event occurs in that grid point. The number of days with daily minimum temperatures below the threshold in the winter of each year (DJF) for each grid point is defined as the extreme low temperature days (ELDs). All grid points define their respective ELDs in the same way.

2.3. Data analysis methods

The statistical methods used in this study include sliding average, cumulative anomaly, empirical orthogonal function decomposition (EOF) [20], Pearson's correlation analysis, and lead-lag regression and correlation. North's test was used for the results of EOF decomposition [21] and Students' t-test was used for correlation and regression analysis [22].

3. Results and discussion

3.1. The spatiotemporal distribution of ELDS

First, this paper studies the spatiotemporal distribution characteristics of winter extreme low temperature days in China. The threshold pattern exhibits a pronounced "south-high, north-low" gradient (Fig. 1a), with the positive and negative threshold boundary line located at 25°N.[8,9]. ELDs demonstrate nationwide uniformity, averaging 17-18 days annually (Fig. 1b). Notably, the annual mean daily minimum temperature increased by about 3°C, rising from -14°C in the 1960s to -11°C in recent decades. Especially after 1985, events with temperatures below -13.5°C have significantly decreased, and events below -13°C have nearly disappeared after 2000. Concurrently, ELDs decreased by over 50%, declining from 10-12 days/year (1960s) to 4-6 days/year (2010s), with heightened interannual variability during 1960s-1980s followed by stabilization post-1985. Cumulative anomaly analysis reveals that there was an interdecadal abrupt change in ELDs in 1985 [8,10,23] (Fig. 1c). The annual average minimum temperature rise and the decrease in ELDs occur synchronously, showing a significant negative correlation feature.

Figure 1. (a) Spatial distribution of the 5th-percentile threshold for extreme low temperatures (°C) in China; (b) Distribution of winter extreme low-temperature days (ELDs); (c) Time series (1961-2020) of: annual mean daily minimum temperature (red solid line), its 9-year sliding average (dashed line), ELDs (bars), and cumulative ELD anomaly (blue solid line), in 1961-2020

This study further conducted a phased statistical analysis of ELDs in December, January, and February in China every 20 years since 1961, and calculated the annual time series of the number of ELDs for each month (Fig. 2). Horizontal comparisons reveals that there is a monthly variation in the high-value areas of ELDs, with a concentration in December in northern, southern and southwestern China; a winter peak was reached in January; and shifted to the northwestern part in February (Figs. 2a-i). Vertical comparisons show that the spatial distribution of ELDs within the same month exhibits intergenerational stability, but the number of ELDs “gradually decreases” in the later period compared to the earlier period, which confirms that the increase in winter average temperature under the background of global warming may have an inhibiting effect on extreme low temperature events, and that regional response differences are significant (with greater attenuation in the northeast and the Qinghai-Tibet Plateau region) (Fig. 2a-i).Meanwhile, except for a few years, the average ELDs in each winter month has significantly decreased since the mid-1980s (Fig. 2j). The above spatiotemporal features corroborate the abrupt climate change phenomenon of ELDs around the 1980s [10,12,19], and the warming trend is more pronounced in the north.

This study identifies the principal modes of extreme low-temperature days (ELDs) through EOF analysis [20] The leading mode (EOF1) exhibits a nationwide coherent pattern with maximum loadings over central China (Fig. 3a), consistent with climatological ELDs distribution (Fig. 1b). EOF2 displays a north-south dipole structure where ELDs in northwestern/northeastern/northern China inversely correlate with southern/southwestern regions (Fig. 3c). These main spatial modes are consistent with previous studies [24-26]. EOF1 and EOF2 together explain 62.52% of the total variance, which can well explain the spatiotemporal characteristics of ELDs. Cumulatively explaining 62.52% variance, both modes meet North et al. [27] criteria for statistical independence, indicating distinct physical drivers.

The principal component of EOF1 (PC1) reveals an interdecadal regime shift transitioning from negative to positive phase around 1985 (Fig. 3b), consistent with China's ELD transition (Fig. 1c) [10]. EOF2's component (PC2) shows prominent interannual variability during 1961–2023 (Fig. 3d), with amplitude reduction post-1980s.

3.2. The correlation between ECDs and atmospheric circulation

To elucidate the physical drivers of these modes, contemporaneous correlations are made between the principal components and the global sea level pressure, the 850 hPa wind field, and the 500 hPa versus 200 hPa geopotential height field and wind field respectively (Fig. 4).

For EOF1’s regionally consistent negative pattern, coherent negative geopotential height anomalies over Arctic/northern Siberia at 200–500 hPa accompany cyclonic circulation anomalies (Fig. 4a-b), suggesting polar vortex expansion. This configuration facilitates southward cold air intrusion and amplifies meridional perturbations in the westerly jet stream [28]. Concurrently, mid-low latitudes exhibit widespread positive height anomalies (Fig. 4a-b). At surface levels, Siberia shows depressed sea-level pressure with 850 hPa cyclonic anomalies, while the Northwest Pacific displays elevated pressure and anticyclonic circulation (Fig. 4c).

The atmospheric circulation associated with the north-south dipole distribution (Fig. 3c) indicates that the winter-time geopotential height and wind field exhibit a quasi-baric distribution from bottom to top (Figures 4d–f). At upper-Mid Troposphere (200–500 hPa), strong positive height anomalies over the Arctic–West Siberia juxtaposed with negative anomalies across northern China–Mongolian Plateau, accompanied by cyclonic circulation (Fig. 4d-e). This configuration reflects polar vortex intensification/southward displacement, enhancing meridional jet stream perturbations that channel cold air into mid-latitudes – a signature resembling the negative Arctic Oscillation phase (AO-) [29,30]. Ural blocking high development promotes Siberian cold air accumulation, with cyclonic steering guiding westerly cold surges into northern China (Fig. 4e). Simultaneously, 850 hPa cyclonic anomalies generate southerly flow over southern China, inhibiting cold air penetration southward (Fig. 4f). The transition from low-level (850 hPa) convergence to upper-level (200 hPa) divergence creates a coupled vertical profile conducive to convective intensification and extreme cold/precipitation events.

3.3. Contemporaneous correlation between ECDs and SSTA

Aiming to further investigate the physical connection between sea surface temperature anomaly（SSTA）and ELDs, the contemporaneous correlation distributions of PC1 and PC2 with global SSTA are calculated (Fig. 5).

For EOF1, the North Pacific Ocean shows a positive-negative-positive correlation distribution from north to south, with a significant negative anomaly signal near the Aleutian Islands, the largest positive anomaly region occurs near the East China Sea along the West Pacific coast, which exhibits an inverse-phase oscillation throughout the year with the negative anomaly in the central part of the region. And almost all the rest of the ocean area exhibits significant positive anomaly correlation except for the Atlantic Ocean south of Greenland, where a negative anomaly center is observed (Fig. 5a). This large-scale positive correlations may reflect the suppression of extreme low temperature events by the overall trend of SST warming under the background of global warming [31]. For the distribution of the contemporaneous correlation between EOF2 and SSTA (Fig. 5b), a significant positive correlation is observed along the southeast coast of China (0°-30°N, 120°E-180°W), and a La Niña-like SST pattern [32], whereas the correlation with other sea areas is not significant.

3.4. Lead-lag regression of ECDs and SSTA

To further explore the signal of SSTA on the ELDs, this paper calculates the lead lag regressions of PC1 and PC2 with SSTA for each season (Fig. 6). For EOF1(nationally coherent mode), weak negative anomalies emerge in the central North Pacific during preceding summer, flanked by positive anomalies in Alaska Bay and the San Francisco coast (Fig. 6a). These signals intensify seasonally: Alaskan positive anomalies expand in fall while central Pacific negative anomalies weaken (Fig. 6b), culminating in amplified phase-locked patterns by late spring (Fig. 6d). Notably, Indian and Atlantic Ocean SSTAs exhibit persistent negative correlations with ELDs across all seasons [31].

EOF2(north-south dipole mode) reveals contrasting phase evolution: its North Pacific precursor anomalies oppose EOF1 during pre-summer to fall (Figs. 6e-f vs. 6a-b) but align during winter to spring (Figs. 6g-h vs. 6c-d). It is also noted that the signal of "east-negative-west-positive" dipole variations in the equatorial Pacific (5°S-5°N) gradually strengthens from the pre-summer to the later spring.

4. Conclusion

This study investigates winter extreme low-temperature events in China, analyzing the spatiotemporal distribution of extreme low-temperature days (ELDs) and their physical linkages with atmospheric circulation and sea surface temperature anomalies (SSTA). Key findings are:

(1) Spatiotemporal Evolution of ELDs: From 1961 to 2020, winter ELDs in China declined by over 50%, with minimum temperatures rising from −14°C to −11°C and sub−13°C extremes disappearing post-2000; EOF analysis reveals two dominant modes explaining 62.52% of the variance: a nationwide decline mode (EOF1, 48.63%) with an interdecadal shift around 1985, and a north–south dipole mode (EOF2, 13.89%) with weakened interannual fluctuations post-1980s.

(2) Atmospheric Circulation Linkages: EOF1 is linked to negative geopotential height anomalies and southward-expanding polar vortices that enhance cold air intrusion but are offset by mid–low latitude blocking highs, while EOF2 features a quasi-barotropic AO-negative-like structure with Arctic–West Siberian highs and northern China lows driving cold surges into northern China and creating a thermal dipole via vertical circulation coupling.

(3) SSTA Connections: EOF1 shows broad positive correlations with winter SSTA across most oceans and year-round tropical (10°S–10°N) SSTs, while EOF2 is associated with La Niña–like tropical Pacific dipole SSTs that intensify from pre-summer to post-spring, though their tropical linkages remain seasonally variable and statistically insignificant.

Limitations and Future Work: Results are limited to correlative analyses; physical processes and causality require numerical experiments. External forcings like sea ice, Eurasian snowpack [33,34], and other factors [35], which modulate events via surface albedo/energy balance [5] remain undiscussed. Future studies should employ high-resolution regional climate models to explore multi-circle interactions, improving predictions of extreme cold events under climate change.